# Data manipulation packagelibrary(tidyverse) # Easily Install and Load the 'Tidyverse'library(psych) # Procedures for Psychological, Psychometric, and Personality Research# Modelling packageslibrary(nlme) # Linear and Nonlinear Mixed Effects Modelslibrary(betareg) # Beta Regressionlibrary(glmmTMB) # Generalized Linear Mixed Models using Template Model Builderlibrary(AICcmodavg) # Model Selection and Multimodel Inference Based on (Q)AIC(c)library(MuMIn) # Multi-Model Inference#Modelling utility packageslibrary(broom) # Convert Statistical Objects into Tidy Tibbleslibrary(broom.mixed) # Tidying Methods for Mixed Modelslibrary(insight) # Easy Access to Model Information for Various Model Objectslibrary(modelbased) # Estimation of Model-Based Predictions, Contrasts and Meanslibrary(parameters) # Processing of Model Parameterslibrary(performance) # Assessment of Regression Models Performancelibrary(DHARMa) # Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Modelslibrary(ggeffects) # Create Tidy Data Frames of Marginal Effects for 'ggplot' from Model Outputs# Grapich packageslibrary(GGally) # Extension to 'ggplot2'library(patchwork) # The Composer of Plots# Table output packageslibrary(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntaxlibrary(DT)library(showtext)font_add_google("Roboto", "Roboto")showtext_auto()showtext_opts(dpi =300)
Variables description
Table 1: List of variables used to modelling spatial co-occurrence patterns of competitive interactions
Variable
Description
Units
Source
SIF
Spatial species interaction factor. SIF<1: segregation, SIF=1: independence, SIF>1= aggregation
Average distance (km) between camera trap stations reported. Only used for spatial co-occurrence models
Average liner distance in km
Each study
Samp_dur
duration of sampling in months. Only used for temporal co-occurrence models
Months
Each study
Locality
Each sampling location where species co-occurrence information was obtained, within each study
category
Each study
Label
Each study unique ID
category
Each study
Data
The data correspond to studies evaluating the spatial and temporal co-occurrence of mammals of the order Carnivora. In all cases the species are of the same order and present ecological interactions of competition.
Spatial co-occurrence data of carnivorous mammals were measured using the species interaction factor (SIF). This is a parameter derived from multi-species occupancy models (Waddle et al. 2010; Richmond, Hines, and Beissinger 2010).
Code
# Load the spatial co-ocurrence data base of competitive interactionsspatC_db <-read_delim("Data/Model data/spatC_db.csv", delim =";" )[,-1] %>%# Select the variablesselect(SIF, D_competitor, D_family, S_competitor, S_family, mass_ratio, diet_dist, Lat, Lat_abs, p_distance, Avg_dist, Locality, Label)DT::datatable(spatC_db)
The Table 2 contain the general description of the numeric variables of spatial data base. The table was constructed with psych package (Revelle 2022)
Code
Spatnum_summary <- spatC_db %>%select_if(is.numeric) %>%describe(. ,fast = T) kbl(Spatnum_summary, caption ="General description of the spatial co-occurrence data base for competitive interactions", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)
Table 2: General description of the spatial co-occurrence data base for competitive interactions.
vars
n
mean
sd
min
max
range
se
SIF
1
216
1.028
0.843
-6.195
5.610
11.805
0.057
mass_ratio
2
211
1.260
1.012
0.027
4.817
4.790
0.070
Lat
3
216
12.484
30.095
-41.234
58.513
99.747
2.048
Lat_abs
4
216
30.170
12.158
0.553
58.513
57.960
0.827
p_distance
5
166
0.145
0.045
0.037
0.200
0.163
0.003
Avg_dist
6
191
1.699
1.418
0.250
10.600
10.350
0.103
We use the GGAlly package to explore the visual relationship between numerical variables (Schloerke et al. 2021). We used pearson’s correlation coefficient
The temporal data corresponds to the co-occurrence of Carnivora mammals measured by overlap coefficient of kernel activity curves (Ridout and Linkie 2009).
The Table 3 contain the general description of the numeric variables of spatial data base.
Code
Tempnum_summary <- tempC_db %>%select_if(is.numeric) %>%describe(. ,fast = T) kbl(Tempnum_summary, caption ="General description of the temporal co-occurrence data base for competitive interactions", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)
Table 3: General description of the temporal co-occurrence data base for competitive interactions.
vars
n
mean
sd
min
max
range
se
Ov_coeff
1
1274
0.636
0.215
0.060
0.970
0.910
0.006
mass_ratio
2
1267
1.408
1.026
0.000
5.678
5.678
0.029
Lat
3
1274
12.429
25.640
-41.234
58.539
99.773
0.718
Lat_abs
4
1274
24.814
13.992
0.553
58.539
57.986
0.392
p_distance
5
915
0.149
0.044
0.037
0.212
0.175
0.001
Samp_dur
6
1015
19.102
29.089
0.780
146.000
145.220
0.913
We use the GGAlly package to explore the visual relationship between numerical variables (Schloerke et al. 2021). We used pearson’s correlation coefficient
The numerical covariates were standardized to mean 0 and standard deviation 1, which facilitates model fitting and direct comparison of regression coefficients. In addition, pairs of species that do not contain information on any of the variables (cells with NAs) were eliminated.
Because we aimed to evaluate whether the patterns found varied depending on the dominant competitor family, we made subsets of the databases. For spatial co-occurrence we only subset the felidae family. For temporal co-occurrence we subset the families Felidae, Canidae and Mustelidae.
Code
# Function to standardized the numeric variablesscale_this <-function(data){scale(data) %>%as.vector() }# Select the variables for Spatial and temporal data setsSC_vars <-c("mass_ratio", "Lat_abs", "p_distance", "Avg_dist")TC_vars <-c("mass_ratio", "Lat_abs", "p_distance", "Samp_dur")# All data SpatialspatC_modsdf <- spatC_db %>%drop_na() %>%mutate(across(all_of(SC_vars), scale_this)) dim(spatC_modsdf)
[1] 149 13
Code
#Subset for spatial data# Felidae SpatialspatC_modsdf_F <- spatC_modsdf %>%filter(D_family =="Felidae")dim(spatC_modsdf_F)
[1] 85 13
Code
# All data TemporaltempC_modsdf <- tempC_db %>%drop_na() %>%mutate(across(all_of(TC_vars), scale_this))dim(tempC_modsdf)
[1] 726 13
Code
#subset for temporal data# Felidae TemporaltempC_modsdf_F <- tempC_modsdf %>%filter(D_family =="Felidae")dim(tempC_modsdf_F)
Identify the presence of outliers or observations that may influence the models. For the spatial co-occurrence data we fit a general linear model with all covariates and interactions. We then checked for the presence of extreme data using Laverage distance plots from performance package (Lüdecke et al. 2021).
Code
# Fit model with all dataSCmods_outl <-lm(SIF ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2, data= spatC_modsdf)# Fit model with only felids as dominantSCmods_outl_F <-lm(SIF ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2, data= spatC_modsdf_F)# Check outliers for each modelSC_out <-check_outliers(SCmods_outl)SC_out_F <-check_outliers(SCmods_outl_F)# Leverage plot for all data modelplot(SC_out)+theme_bw()+theme(text=element_text(family ="Roboto"))# Leverage plot for felidae as dominant competitorplot(SC_out_F)+theme_bw()+theme(text=element_text(family ="Roboto"))
Outliers of all data spatial model
Outliers of felidae dominant competior spatial model
Spatial co-occurrence leverage distance
Collinearity
We checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value > 10 indicates a high collinearity of the variables (Alain F. Zuur, Ieno, and Elphick 2010).
#| message: false#| warning: false#| label: fig-Scoll#| fig-cap: Spatial co-occurrence collinearity#| fig-subcap: #| - "VIF of all data spatial model"#| - "VIF of felidae dominant competior spatial model"#| layout-ncol: 2#| fig-width: 8# Fit model with all dataspatC_coll <-lm(SIF ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist, data= spatC_modsdf)# Fit model with only felids as dominantspatC_coll_F <-lm(SIF ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist, data= spatC_modsdf_F)# Check outliers for each modelspatC_colltab <-check_collinearity(spatC_coll)spatC_colltab_F <-check_collinearity(spatC_coll_F)# Vif plotsplot(spatC_colltab)+labs(subtitle ="All data Spatial model")+theme_bw()+theme(legend.position ="none",text=element_text(family ="Roboto"))
We first check the assumptions of the model. To do this we fit the most complex model containing all available variable interactions and random factors. We then visually evaluated the normality of the residuals, linearity and homogeneity of variance using the performance package (Lüdecke et al. 2021).
We detect deviations from assumptions. Particularly there is a kurtosis in the distribution of the residuals. There are also deviations from linearity. For that reason, we fit a T-student family glm using glmmTMB package (Brooks et al. 2017). To evaluate the GLMM we use the DHARMa package (Hartig 2021), which applies simulated residuals to evaluate the assumptions.
$uniformity
$uniformity$details
catPred: Diff_diet
One-sample Kolmogorov-Smirnov test
data: dd[x, ]
D = 0.23074, p-value = 0.2802
alternative hypothesis: two-sided
------------------------------------------------------------
catPred: Same_diet
One-sample Kolmogorov-Smirnov test
data: dd[x, ]
D = 0.10399, p-value = 0.4252
alternative hypothesis: two-sided
$uniformity$p.value
[1] 0.2802263 0.4251758
$uniformity$p.value.cor
[1] 0.5604526 0.5604526
$homogeneity
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.3309 0.5667
83
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.01022
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.09036
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.07346
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.09331
alternative hypothesis: both
Uniformity of variance test of Diet category
Uniformity of variance test of mass ratio
Uniformity of variance test of absolute latitude
Uniformity of variance test of phylogenetic distance
Uniformity of variance test of Average distance
We detected non-uniformity in the range of body mass ratio, phylogenetic distance, absolute latitude and average camera distance. Additionally, we identified that observation 84 (previously identified as outlier) does have an effect on model fit. To improve the fit, we modeled the variation of the variables by means of the dispformula term, from the glmm TMB package (Brooks et al. 2017).
Following the Alain F. Zuur et al. (2009a) protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the previously selected model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package (Mazerolle 2023).
Code
#| tbl-cap: Random-effects structure selection table for spatial co-occurrence data# Fit model without random effectsSpat_r0 <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2,data=spatC_modsdf, family=t_family(link ="identity"), REML = T)# Fit model with label as random effectSpat_r1 <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label),data=spatC_modsdf, family=t_family(link ="identity"), REML = T)# Fit model with label and Locality as random effectsSpat_r2 <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label/Locality),data=spatC_modsdf, family=t_family(link ="identity"), REML = T)# Rank the models with the AICspat_AICtab <-aictab(cand.set =list(Spat_r0, Spat_r1, Spat_r2),modnames =c("no random","Label random", "Label/Locality random"),second.ord = F)kbl(spat_AICtab, caption ="Random-effects structure selection table for spatial co-occurrence data", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)
Random-effects structure selection table for spatial co-occurrence data
Modnames
K
AIC
Delta_AIC
ModelLik
AICWt
LL
Cum.Wt
2
Label random
25
293.925
0.000
1.000
0.496
-121.962
0.496
1
no random
24
294.786
0.861
0.650
0.322
-123.393
0.818
3
Label/Locality random
26
295.925
2.000
0.368
0.182
-121.962
1.000
The ?@tbl-Sraic suggests that both models with Label as random effects and without random effects are equally supported. To visualize the variation of the random parameters of each group, we use the estimate_grouplevel function of the modelbased package (Makowski et al. 2020).
Code
estimate_grouplevel(Spat_r1) %>%# Get random parametersplot()+#plot# aesthetic changestheme_bw()+theme(text=element_text(family ="Roboto"))
There is no significant variation from each group level in random effects. We now check the goodness of fit of these models.
Random-effects structure selection table for spatial co-occurrence data when Felidae is dominant competitor
Modnames
K
AIC
Delta_AIC
ModelLik
AICWt
LL
Cum.Wt
no random
27
180.148
NA
NA
NA
-63.074
NA
Label random
28
NA
NA
NA
NA
NA
NA
Label/Locality random
29
183.280
NA
NA
NA
-62.640
NA
In the case of Felidae the model selection indicates that random effects are not adequate. This is because there is no variation between Label or Locality levels.
All label coefficient are the same for felids data
Fixed predictor selections
To select the fixed effects structure we also used the Akaike information criterion (AIC). Models with a delta AIC <2 are considered equally plausible (Burnham and Anderson 2002). Because we are going to choose the fixed effects structure, we refit the previously identified model without using the restricted maximum likelihood (REML) (Alain F. Zuur et al. 2009b). Since co-occurrence patterns may be generated by the interaction of variables, we fit all possible combinations of variables, limiting to a maximum of three parameters per model. To do this we will used the dredge function of the MuMIn package (Barton 2020).
All data
Code
# Create a global modelSC_global <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2,data=spatC_modsdf, family=t_family(link ="identity"), REML = F,na.action ="na.fail",start =list(psi =log(1.69)),map =list(psi =factor(NA)))#Fit all posible combination limit to 3 parametersSC_selec <-dredge(SC_global, rank ="AIC", m.lim=c(NA,3))
Fixed-effects structure selection table for spatial co-occurrence data
cond((Int))
disp((Int))
cond(Avg_dist)
cond(diet_dist)
cond(Lat_abs)
cond(mass_ratio)
cond(I(mass_ratio^2))
cond(p_distance)
cond(Avg_dist:diet_dist)
cond(Avg_dist:Lat_abs)
cond(Avg_dist:mass_ratio)
cond(Avg_dist:I(mass_ratio^2))
cond(Avg_dist:p_distance)
cond(diet_dist:Lat_abs)
cond(diet_dist:mass_ratio)
cond(diet_dist:I(mass_ratio^2))
cond(diet_dist:p_distance)
cond(Lat_abs:mass_ratio)
cond(Lat_abs:p_distance)
cond(mass_ratio:p_distance)
cond(I(mass_ratio^2):Lat_abs)
cond(I(mass_ratio^2):mass_ratio)
cond(I(mass_ratio^2):p_distance)
df
logLik
AIC
delta
weight
49
1.032
+
NA
NA
NA
NA
-0.024
0.028
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-75.509
159.017
0.000
0.095233872
1048625
1.029
+
NA
NA
NA
NA
-0.016
0.039
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.018
5
-74.654
159.308
0.291
0.082332924
17
1.026
+
NA
NA
NA
NA
-0.021
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
-76.864
159.727
0.710
0.066762335
50
1.033
+
-0.009
NA
NA
NA
-0.024
0.026
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-75.370
160.741
1.724
0.040225789
51
1.025
+
NA
+
NA
NA
-0.023
0.030
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-75.476
160.952
1.935
0.036186594
57
1.030
+
NA
NA
NA
-0.006
-0.021
0.028
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-75.481
160.962
1.945
0.036005284
18
1.028
+
-0.015
NA
NA
NA
-0.022
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-76.496
160.992
1.975
0.035481985
53
1.032
+
NA
NA
-0.002
NA
-0.024
0.028
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-75.502
161.005
1.988
0.035249935
41
1.010
+
NA
NA
NA
-0.030
NA
0.026
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-76.551
161.102
2.085
0.033578960
19
1.042
+
NA
+
NA
NA
-0.022
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-76.682
161.364
2.347
0.029451872
9
1.007
+
NA
NA
NA
-0.026
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
-77.728
161.457
2.440
0.028121123
131113
1.013
+
NA
NA
NA
-0.027
NA
0.023
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.026
NA
NA
NA
5
-75.796
161.591
2.574
0.026291297
1
1.005
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
2
-78.829
161.659
2.642
0.025415233
25
1.025
+
NA
NA
NA
-0.004
-0.019
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-76.851
161.703
2.686
0.024868031
21
1.026
+
NA
NA
0.002
NA
-0.021
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-76.856
161.712
2.695
0.024752479
33
1.008
+
NA
NA
NA
NA
NA
0.023
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
-77.922
161.845
2.828
0.023161029
43
0.987
+
NA
+
NA
-0.031
NA
0.035
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.154
162.308
3.291
0.018368064
8211
1.030
+
NA
+
NA
NA
-0.014
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
5
-76.207
162.414
3.397
0.017427872
10
1.007
+
-0.018
NA
NA
-0.031
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-77.253
162.506
3.489
0.016639325
20
1.044
+
-0.015
+
NA
NA
-0.024
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.306
162.613
3.596
0.015775034
524313
1.037
+
NA
NA
NA
-0.026
-0.048
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.015
NA
5
-76.324
162.648
3.631
0.015501177
42
1.010
+
-0.013
NA
NA
-0.033
NA
0.024
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.327
162.654
3.637
0.015453356
530
1.030
+
-0.012
NA
NA
NA
-0.027
NA
NA
NA
NA
-0.008
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.416
162.832
3.815
0.014139679
26
1.025
+
-0.017
NA
NA
-0.009
-0.019
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.427
162.854
3.837
0.013986375
45
1.010
+
NA
NA
-0.007
-0.031
NA
0.028
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.472
162.943
3.926
0.013371798
22
1.028
+
-0.015
NA
0.003
NA
-0.022
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.478
162.956
3.939
0.013285763
23
1.046
+
NA
+
-0.005
NA
-0.023
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.655
163.310
4.293
0.011133483
262165
1.026
+
NA
NA
-0.003
NA
-0.020
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.006
NA
NA
5
-76.670
163.341
4.323
0.010963681
35
0.990
+
NA
+
NA
NA
NA
0.030
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-77.674
163.348
4.331
0.010921479
27
1.042
+
NA
+
NA
0.001
-0.022
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.682
163.364
4.347
0.010837038
2
1.005
+
-0.009
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
-78.697
163.394
4.377
0.010671899
13
1.006
+
NA
NA
-0.002
-0.027
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-77.722
163.445
4.428
0.010407494
11
1.008
+
NA
+
NA
-0.026
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-77.726
163.452
4.435
0.010369519
5
1.005
+
NA
NA
0.003
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
-78.818
163.635
4.618
0.009461142
3
1.008
+
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
-78.823
163.645
4.628
0.009413222
29
1.025
+
NA
NA
0.002
-0.003
-0.019
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-76.847
163.694
4.677
0.009185268
34
1.008
+
-0.004
NA
NA
NA
NA
0.022
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-77.897
163.795
4.778
0.008735579
37
1.008
+
NA
NA
-0.001
NA
NA
0.023
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-77.922
163.843
4.826
0.008526966
16419
0.962
+
NA
+
NA
NA
NA
0.078
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
5
-76.957
163.914
4.897
0.008228560
65573
1.005
+
NA
NA
0.008
NA
NA
0.027
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.024
NA
NA
NA
NA
5
-77.021
164.041
5.024
0.007722202
14
1.007
+
-0.018
NA
-0.002
-0.032
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.248
164.496
5.479
0.006151739
266
1.007
+
-0.018
NA
NA
-0.031
NA
NA
NA
NA
-0.001
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.252
164.504
5.487
0.006126574
12
1.008
+
-0.018
+
NA
-0.031
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.253
164.506
5.489
0.006123161
4107
1.009
+
NA
+
NA
-0.039
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.421
164.842
5.825
0.005175900
39
0.985
+
NA
+
0.008
NA
NA
0.031
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.602
165.204
6.187
0.004319037
32781
1.008
+
NA
NA
-0.003
-0.026
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.007
NA
NA
NA
NA
NA
5
-77.646
165.291
6.274
0.004134075
36
0.990
+
-0.003
+
NA
NA
NA
0.029
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.662
165.323
6.306
0.004067907
6
1.005
+
-0.009
NA
0.004
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-78.678
165.356
6.339
0.004002483
4
1.007
+
-0.009
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-78.692
165.384
6.367
0.003946027
15
1.010
+
NA
+
-0.004
-0.027
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.711
165.423
6.406
0.003870416
134
1.007
+
0.001
NA
0.011
NA
NA
NA
NA
-0.017
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.766
165.532
6.515
0.003664890
7
1.006
+
NA
+
0.002
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
-78.817
165.633
6.616
0.003484069
1058
1.009
+
-0.003
NA
NA
NA
NA
0.022
NA
NA
NA
NA
0.002
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.894
165.788
6.771
0.003225380
38
1.008
+
-0.004
NA
0.000
NA
NA
0.022
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-77.897
165.795
6.778
0.003213936
2055
0.985
+
NA
+
0.041
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-78.513
167.026
8.009
0.001736160
68
1.007
+
-0.016
+
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-78.569
167.139
8.122
0.001641099
8
1.006
+
-0.009
+
0.004
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
-78.678
167.356
8.339
0.001472433
Felidae
Note that this time the limit number of parameters in the function is 7. This is due to the number of parameters of the dispersion terms, which add 4 additional parameters.
Code
spatC_modsdf_F <- spatC_modsdf_F[-84,]SC_global_F <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2,dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F, family=t_family(link ="identity"), REML = F,start =list(psi =log(2.87)),map =list(psi =factor(NA)),na.action ="na.fail")SC_selec_F <-dredge(SC_global_F, rank ="AIC", fixed =c("disp(Avg_dist)", "disp(Lat_abs)", "disp(mass_ratio)", "disp(p_distance)"), # Mantain dispersion parameters in all modelsm.lim=c(NA,7)) # 7 parameteres becasue we must add the dispersion terms
Code
kbl(SC_selec_F, caption ="Fixed-effects structure selection table for spatial co-occurrence data when Felidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)
Fixed-effects structure selection table for spatial co-occurrence data when Felidae is dominant competitor
cond((Int))
disp((Int))
cond(Avg_dist)
cond(diet_dist)
cond(Lat_abs)
cond(mass_ratio)
cond(I(mass_ratio^2))
cond(p_distance)
cond(Avg_dist:diet_dist)
cond(Avg_dist:Lat_abs)
cond(Avg_dist:mass_ratio)
cond(Avg_dist:I(mass_ratio^2))
cond(Avg_dist:p_distance)
cond(diet_dist:Lat_abs)
cond(diet_dist:mass_ratio)
cond(diet_dist:I(mass_ratio^2))
cond(diet_dist:p_distance)
cond(Lat_abs:mass_ratio)
cond(Lat_abs:p_distance)
cond(mass_ratio:p_distance)
cond(I(mass_ratio^2):Lat_abs)
cond(I(mass_ratio^2):mass_ratio)
cond(I(mass_ratio^2):p_distance)
disp(Avg_dist)
disp(Lat_abs)
disp(mass_ratio)
disp(p_distance)
df
logLik
AIC
delta
weight
262165
1.041
-1.609
NA
NA
-0.001
NA
-0.107
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.067
NA
NA
-0.033
0.109
0.024
-0.183
9
-29.918
77.836
0.000
0.106482503
21
0.995
-1.574
NA
NA
-0.041
NA
-0.034
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.036
0.147
0.025
-0.138
8
-31.103
78.206
0.370
0.088479291
17
1.018
-1.581
NA
NA
NA
NA
-0.030
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.038
0.127
0.018
-0.086
7
-32.233
78.466
0.630
0.077718447
1048625
1.011
-1.613
NA
NA
NA
NA
-0.016
0.024
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.024
0.096
0.125
-0.053
-0.068
9
-30.878
79.756
1.920
0.040774185
22
0.997
-1.584
-0.019
NA
-0.043
NA
-0.036
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.038
0.146
0.033
-0.151
9
-30.971
79.943
2.107
0.037134205
53
0.995
-1.576
NA
NA
-0.041
NA
-0.034
0.009
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.038
0.159
0.033
-0.146
9
-31.035
80.069
2.233
0.034858667
23
1.011
-1.575
NA
+
-0.046
NA
-0.034
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.039
0.148
0.024
-0.148
9
-31.066
80.132
2.296
0.033781932
19
0.987
-1.577
NA
+
NA
NA
-0.030
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.019
0.135
0.025
-0.092
8
-32.079
80.157
2.321
0.033358046
29
0.995
-1.574
NA
NA
-0.041
0.000
-0.034
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.035
0.147
0.025
-0.138
9
-31.103
80.206
2.370
0.032549878
25
1.023
-1.589
NA
NA
NA
0.018
-0.037
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.017
0.125
0.019
-0.104
8
-32.108
80.216
2.380
0.032398022
524313
1.050
-1.598
NA
NA
NA
-0.009
-0.091
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.023
NA
0.034
0.104
0.042
-0.120
9
-31.185
80.369
2.533
0.030005702
18
1.019
-1.590
-0.012
NA
NA
NA
-0.031
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.052
0.118
0.023
-0.088
8
-32.185
80.370
2.534
0.029989393
49
1.020
-1.584
NA
NA
NA
NA
-0.031
0.007
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.039
0.132
0.023
-0.092
8
-32.189
80.379
2.543
0.029863756
530
1.047
-1.604
0.040
NA
NA
NA
-0.080
NA
NA
NA
NA
-0.064
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.100
0.084
0.028
-0.111
9
-31.241
80.482
2.646
0.028365000
13
0.962
-1.519
NA
NA
-0.049
-0.048
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.052
0.147
0.013
-0.065
8
-32.453
80.906
3.070
0.022942770
1
0.982
-1.543
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.047
0.099
0.015
0.004
6
-34.498
80.997
3.161
0.021923939
9
0.990
-1.537
NA
NA
NA
-0.032
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.102
0.113
-0.003
0.001
7
-33.736
81.472
3.636
0.017286222
32781
0.983
-1.550
NA
NA
-0.044
-0.014
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.044
NA
NA
NA
NA
NA
0.103
0.114
0.015
-0.063
9
-31.806
81.612
3.777
0.016113987
5
0.963
-1.527
NA
NA
-0.032
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.023
0.135
0.039
-0.042
7
-33.864
81.729
3.893
0.015203657
51
0.976
-1.585
NA
+
NA
NA
-0.031
0.014
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.017
0.146
0.038
-0.104
9
-31.922
81.843
4.008
0.014356678
8211
0.967
-1.573
NA
+
NA
NA
-0.016
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
0.027
0.146
-0.006
-0.094
9
-31.922
81.845
4.009
0.014348124
2055
0.880
-1.546
NA
+
0.167
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.041
0.156
0.083
-0.067
9
-31.953
81.907
4.071
0.013909338
27
0.997
-1.584
NA
+
NA
0.013
-0.035
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.010
0.132
0.025
-0.103
9
-32.020
82.040
4.204
0.013010901
20
0.989
-1.585
-0.011
+
NA
NA
-0.031
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.031
0.128
0.029
-0.095
9
-32.040
82.080
4.244
0.012756341
3
0.927
-1.533
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.031
0.120
0.010
0.008
7
-34.074
82.148
4.312
0.012327659
57
1.023
-1.591
NA
NA
NA
0.016
-0.037
0.005
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.021
0.128
0.023
-0.107
9
-32.083
82.167
4.331
0.012213911
14
0.957
-1.533
-0.040
NA
-0.059
-0.064
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.031
0.157
0.032
-0.088
9
-32.091
82.183
4.347
0.012116875
26
1.023
-1.591
-0.004
NA
NA
0.016
-0.037
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.024
0.122
0.021
-0.103
9
-32.104
82.208
4.372
0.011966312
50
1.020
-1.591
-0.011
NA
NA
NA
-0.032
0.006
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.051
0.124
0.027
-0.093
9
-32.154
82.307
4.471
0.011386440
11
0.933
-1.530
NA
+
NA
-0.032
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.070
0.132
0.003
-0.007
8
-33.306
82.611
4.775
0.009780307
45
0.962
-1.521
NA
NA
-0.051
-0.051
NA
0.011
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.047
0.167
0.025
-0.083
9
-32.361
82.723
4.887
0.009248308
15
0.946
-1.519
NA
+
-0.045
-0.047
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.050
0.147
0.012
-0.057
9
-32.415
82.830
4.994
0.008768419
2
0.983
-1.535
0.012
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.037
0.104
0.004
0.004
7
-34.451
82.901
5.065
0.008460080
16419
0.538
-1.580
NA
+
NA
NA
NA
0.467
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
0.121
0.128
0.080
-0.021
9
-32.469
82.939
5.103
0.008302923
33
0.982
-1.542
NA
NA
NA
NA
NA
-0.002
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.048
0.096
0.012
0.008
7
-34.494
82.987
5.151
0.008103121
266
1.008
-1.589
0.007
NA
NA
-0.007
NA
NA
NA
NA
0.063
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.113
0.101
0.032
-0.072
9
-32.588
83.177
5.341
0.007370482
10
0.990
-1.545
-0.015
NA
NA
-0.036
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.115
0.107
0.004
0.003
8
-33.681
83.362
5.527
0.006717325
41
0.992
-1.538
NA
NA
NA
-0.032
NA
0.004
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.100
0.118
0.000
-0.003
8
-33.726
83.451
5.615
0.006426230
7
0.933
-1.525
NA
+
-0.025
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.011
0.136
0.026
-0.021
8
-33.739
83.478
5.642
0.006340809
131113
1.002
-1.543
NA
NA
NA
-0.037
NA
0.004
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.042
NA
NA
NA
0.197
0.090
-0.084
0.056
9
-32.746
83.492
5.657
0.006294670
6
0.963
-1.523
0.008
NA
-0.031
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.018
0.134
0.030
-0.038
8
-33.843
83.686
5.850
0.005713972
37
0.963
-1.527
NA
NA
-0.032
NA
NA
0.001
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.024
0.136
0.041
-0.044
8
-33.864
83.727
5.891
0.005597270
35
0.919
-1.538
NA
+
NA
NA
NA
0.010
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.022
0.131
0.020
-0.007
8
-34.004
84.008
6.172
0.004864229
4
0.926
-1.527
0.012
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.025
0.124
0.002
0.007
8
-34.027
84.055
6.219
0.004751913
43
0.917
-1.537
NA
+
NA
-0.036
NA
0.019
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.050
0.156
0.022
-0.031
9
-33.076
84.152
6.316
0.004526128
4107
0.932
-1.532
NA
+
NA
-0.053
NA
NA
NA
NA
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
0.064
0.127
0.015
-0.001
9
-33.190
84.380
6.544
0.004039537
12
0.933
-1.538
-0.016
+
NA
-0.037
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.078
0.127
0.008
-0.006
9
-33.240
84.480
6.644
0.003842679
134
0.974
-1.576
-0.024
NA
-0.028
NA
NA
NA
NA
-0.071
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.266
0.074
0.049
-0.047
9
-33.327
84.653
6.817
0.003522759
34
0.982
-1.535
0.012
NA
NA
NA
NA
-0.001
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.038
0.103
0.003
0.005
8
-34.450
84.899
7.063
0.003115193
42
0.991
-1.545
-0.014
NA
NA
-0.036
NA
0.002
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.113
0.110
0.005
0.000
9
-33.677
85.353
7.517
0.002482502
39
0.927
-1.527
NA
+
-0.024
NA
NA
0.007
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.015
0.146
0.033
-0.031
9
-33.703
85.406
7.570
0.002417765
8
0.932
-1.521
0.009
+
-0.024
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.006
0.136
0.017
-0.018
9
-33.709
85.417
7.581
0.002404586
65573
0.961
-1.528
NA
NA
-0.036
NA
NA
-0.004
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.009
NA
NA
NA
NA
-0.048
0.145
0.042
-0.047
9
-33.809
85.617
7.781
0.002175595
38
0.963
-1.524
0.008
NA
-0.031
NA
NA
0.001
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.020
0.136
0.032
-0.040
9
-33.841
85.683
7.847
0.002105420
36
0.918
-1.531
0.014
+
NA
NA
NA
0.012
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.015
0.138
0.013
-0.009
9
-33.935
85.870
8.034
0.001917204
68
0.926
-1.521
0.038
+
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.007
0.130
0.000
0.009
9
-33.992
85.985
8.149
0.001810560
1058
0.984
-1.520
0.023
NA
NA
NA
NA
-0.002
NA
NA
NA
NA
0.021
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.020
0.116
-0.023
0.016
9
-34.341
86.682
8.846
0.001277834
Confidence intervals explorations
Consistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models (Sutherland et al. 2023). Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference (Arnold 2010).
Code
# Function to get 85%ICget_ci <-function(mod_name, mods){ ci_df <- parameters::model_parameters(mods, ci =0.85) %>%select(Parameter, Coefficient, CI_low, CI_high, Component) %>%filter(Component =="conditional") %>%filter(Parameter !="(Intercept)") %>%mutate(Model = mod_name,Informative =case_when( CI_low <0& CI_high <0~"yes", CI_low >0& CI_high >0~"yes",TRUE~"no")) return(ci_df) }# Function to plotci_plot <-function(ci_df){ggplot(ci_df, aes(x=Coefficient, y= Parameter))+geom_pointrange(aes(xmin=CI_low, xmax= CI_high,col= Model, linetype= Informative),position =position_dodge2(0.5), linewidth=1)+scale_linetype_manual(values =c("no"="dashed", "yes"="solid"))+geom_vline(xintercept =0, linetype="dashed")+scale_color_viridis_d()+labs(caption ="estimates with 85% ci intervals",col="Model ID")+theme_bw()+theme(text=element_text(family ="Roboto"))}
Code
# get best models for all data SC_best_mods <-get.models(SC_selec, subset = delta <2)# get best models for felidae subset dataSC_best_mods_F <-get.models(SC_selec_F, subset = delta <2)# Apply the function to obtain the table of coefficients of selected the modelsSC_best_ci <-map2_df(names(SC_best_mods), SC_best_mods, get_ci) SC_best_ci_F <-map2_df(names(SC_best_mods_F), SC_best_mods_F, get_ci) ci_plot(SC_best_ci)+labs(title="All data")ci_plot(SC_best_ci_F)+labs(title="Felidae dominant competitor")
All data best model 85% CI
Felidae data best model 85% CI
Temporal modelling procedure
Influential observations
Identify the presence of outliers or observations that may influence the models. For the temporal co-ocurrence data we fit a beta family generalized linear model with all covariates and interactions (Cribari-Neto and Zeileis 2010). We then checked for the presence of extreme observations using Cook’s distance plots from performance package (Lüdecke et al. 2021).
Outliers of Felidae dominant competior Temporal model
Outliers of Celidae dominant competior Temporal model
Outliers of Mustelidae dominant competior Temporal model
Temporal co-occurrence cook distance distance
A possible outlier was detected for the Mustelidae data as a dominant competitor. It corresponds to observation 27 which exceeded the 0.98 Cook’s distance threshold.
Collinearity
We checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value > 10 indicates a high collinearity of the variables (Alain F. Zuur, Ieno, and Elphick 2010).
Figure 3: Outliers model detection temporal models
Model assumptions
We verify the model assumptions by visual inspection of the simulated residuals from the DHARMa package (Hartig 2021). To do so we fit the more complex model which includes fixed variables and their interactions, as well as random effects.
$uniformity
$uniformity$details
catPred: Diff_diet
One-sample Kolmogorov-Smirnov test
data: dd[x, ]
D = 0.091409, p-value = 0.01375
alternative hypothesis: two-sided
------------------------------------------------------------
catPred: Same_diet
One-sample Kolmogorov-Smirnov test
data: dd[x, ]
D = 0.058897, p-value = 0.1027
alternative hypothesis: two-sided
$uniformity$p.value
[1] 0.01374834 0.10265831
$uniformity$p.value.cor
[1] 0.02749668 0.10265831
$homogeneity
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0 0.9949
724
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.1531
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value < 2.2e-16
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.01388
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value < 2.2e-16
alternative hypothesis: both
Uniformity of variance test of Diet category
Uniformity of variance test of mass ratio
Uniformity of variance test of absolute latitude
Uniformity of variance test of phylogenetic distance
Uniformity of variance test of sampling duration
We detected non-uniformity in the range of diet distance, phylogenetic distance, absolute latitude and mean camera distance. Additionally, we identified that observation 84 (previously identified as outlier) does have an effect on model fit. To improve the fit, we modeled the variation of the variables by means of the dispformula term, from the glmm TMB package (Brooks et al. 2017).
Following the Alain F. Zuur et al. (2009a) protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the previously selected model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package (Mazerolle 2023).
kbl(TC_selec_M, caption ="Fixed-effects structure selection table for temporal co-occurrence data when Mustelidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)
Fixed-effects structure selection table for temporal co-occurrence data when Mustelidae is dominant competitor
(Intercept)
diet_dist
Lat_abs
mass_ratio
I(mass_ratio^2)
p_distance
Samp_dur
diet_dist:Lat_abs
diet_dist:mass_ratio
diet_dist:I(mass_ratio^2)
diet_dist:p_distance
diet_dist:Samp_dur
Lat_abs:mass_ratio
Lat_abs:p_distance
Lat_abs:Samp_dur
mass_ratio:p_distance
mass_ratio:Samp_dur
I(mass_ratio^2):Lat_abs
I(mass_ratio^2):mass_ratio
I(mass_ratio^2):p_distance
I(mass_ratio^2):Samp_dur
p_distance:Samp_dur
df
logLik
AIC
delta
weight
20
-0.051
+
0.387
NA
NA
0.392
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
44.905
-79.811
0.000
7.762372e-01
42
0.550
+
NA
NA
-0.196
NA
-0.368
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
42.156
-74.312
5.499
4.964749e-02
12
0.463
+
0.228
NA
-0.187
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
42.028
-74.056
5.755
4.368093e-02
26
0.510
+
NA
NA
-0.180
0.212
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
41.365
-72.731
7.080
2.251960e-02
50
0.213
+
NA
NA
NA
0.244
-0.367
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
40.879
-71.759
8.052
1.385403e-02
8
0.213
+
0.256
-0.170
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
40.568
-71.135
8.676
1.014190e-02
266
0.690
+
NA
NA
-0.240
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
40.443
-70.886
8.925
8.954648e-03
10
0.651
+
NA
NA
-0.205
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
39.386
-70.772
9.039
8.457681e-03
38
0.310
+
NA
-0.156
NA
NA
-0.367
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
39.992
-69.984
9.827
5.703549e-03
14
0.601
+
NA
-0.090
-0.175
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
39.905
-69.810
10.001
5.228230e-03
27
0.512
NA
0.313
NA
-0.197
0.160
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
39.904
-69.807
10.004
5.221000e-03
11
0.625
NA
0.249
NA
-0.216
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
38.711
-69.423
10.388
4.308490e-03
43
0.619
NA
0.177
NA
-0.215
NA
-0.249
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
39.704
-69.408
10.403
4.276712e-03
36
0.243
+
0.177
NA
NA
NA
-0.257
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
39.603
-69.206
10.605
3.865914e-03
530
0.204
+
NA
NA
NA
0.435
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
39.551
-69.102
10.708
3.670657e-03
4
0.250
+
0.252
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
38.497
-68.994
10.817
3.476237e-03
41
0.739
NA
NA
NA
-0.228
NA
-0.376
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
38.444
-68.889
10.922
3.298852e-03
34
0.341
+
NA
NA
NA
NA
-0.385
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
38.311
-68.622
11.189
2.886762e-03
18
0.297
+
NA
NA
NA
0.258
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
38.178
-68.357
11.454
2.528047e-03
68
0.198
+
0.325
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
39.071
-68.141
11.669
2.270187e-03
15
0.598
NA
0.254
-0.058
-0.198
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
38.934
-67.867
11.943
1.979477e-03
65547
0.589
NA
0.329
NA
-0.178
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.099
NA
NA
NA
NA
5
38.930
-67.860
11.951
1.971698e-03
22
0.303
+
NA
-0.105
NA
0.209
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
38.856
-67.711
12.099
1.830980e-03
134
0.387
+
NA
-0.256
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
38.724
-67.449
12.362
1.605798e-03
524329
0.739
NA
NA
NA
-0.225
NA
-0.451
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.102
NA
5
38.635
-67.270
12.541
1.468064e-03
1058
0.325
+
NA
NA
NA
NA
-0.456
NA
NA
NA
NA
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
38.595
-67.190
12.620
1.411026e-03
57
0.724
NA
NA
NA
-0.225
0.045
-0.374
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
38.552
-67.104
12.707
1.351306e-03
45
0.727
NA
NA
-0.031
-0.219
NA
-0.373
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
38.509
-67.019
12.792
1.295074e-03
4115
0.073
NA
0.540
NA
NA
0.488
NA
NA
NA
NA
NA
NA
NA
-0.322
NA
NA
NA
NA
NA
NA
NA
NA
5
38.386
-66.773
13.038
1.145086e-03
6
0.400
+
NA
-0.171
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
37.341
-66.681
13.129
1.093953e-03
9
0.847
NA
NA
NA
-0.239
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
35.814
-65.628
14.183
6.460790e-04
131085
0.747
NA
NA
0.212
-0.054
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.144
NA
NA
NA
5
37.485
-64.969
14.841
4.647566e-04
2
0.436
+
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
35.422
-64.843
14.967
4.363696e-04
8227
0.458
NA
0.275
NA
NA
NA
-0.350
NA
NA
NA
NA
NA
NA
NA
0.343
NA
NA
NA
NA
NA
NA
NA
5
37.365
-64.730
15.081
4.123065e-04
19
0.281
NA
0.365
NA
NA
0.207
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
36.186
-64.372
15.438
3.448571e-04
25
0.829
NA
NA
NA
-0.235
0.051
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
35.948
-63.896
15.914
2.717865e-04
13
0.830
NA
NA
-0.042
-0.227
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
35.924
-63.847
15.963
2.651965e-04
51
0.293
NA
0.302
NA
NA
0.182
-0.179
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
36.669
-63.338
16.473
2.055938e-04
2055
0.341
NA
0.343
-0.265
NA
NA
NA
NA
NA
NA
NA
NA
0.184
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
36.656
-63.312
16.499
2.029507e-04
23
0.289
NA
0.357
-0.090
NA
0.168
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
36.648
-63.297
16.514
2.013886e-04
7
0.379
NA
0.292
-0.141
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
35.605
-63.209
16.601
1.927815e-04
39
0.375
NA
0.226
-0.131
NA
NA
-0.225
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
36.391
-62.782
17.028
1.557298e-04
3
0.403
NA
0.281
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
34.335
-62.670
17.141
1.472128e-04
35
0.398
NA
0.208
NA
NA
NA
-0.251
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
35.302
-62.605
17.206
1.424715e-04
262169
0.830
NA
NA
NA
-0.227
-0.030
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.087
NA
NA
5
36.106
-62.212
17.599
1.170713e-04
29
0.822
NA
NA
-0.028
-0.228
0.039
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
35.990
-61.981
17.830
1.042925e-04
33
0.525
NA
NA
NA
NA
NA
-0.401
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
33.655
-61.310
18.501
7.458279e-05
37
0.513
NA
NA
-0.115
NA
NA
-0.387
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
34.454
-60.908
18.903
6.099542e-05
32805
0.514
NA
NA
-0.059
NA
NA
-0.361
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.268
NA
NA
NA
NA
NA
5
35.442
-60.884
18.926
6.027970e-05
49
0.506
NA
NA
NA
NA
0.068
-0.397
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
33.883
-59.765
20.046
3.444807e-05
53
0.508
NA
NA
-0.107
NA
0.022
-0.387
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
5
34.473
-58.947
20.864
2.288102e-05
1048625
0.521
NA
NA
NA
NA
-0.001
-0.354
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
-0.21
5
34.280
-58.559
21.252
1.884820e-05
1
0.628
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
2
30.844
-57.687
22.123
1.218964e-05
5
0.612
NA
NA
-0.131
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
31.842
-57.684
22.127
1.216686e-05
17
0.606
NA
NA
NA
NA
0.077
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
3
31.119
-56.239
23.572
5.907693e-06
21
0.606
NA
NA
-0.123
NA
0.024
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA
4
31.864
-55.729
24.082
4.577766e-06
16405
0.611
NA
NA
-0.120
NA
0.048
NA
NA
NA
NA
NA
NA
NA
NA
NA
0.036
NA
NA
NA
NA
NA
NA
5
31.887
-53.774
26.037
1.722455e-06
Confidence intervals explorations
Code
# Function to get 85%ICget_cibeta <-function(mod_name, mods){ ci_df <- parameters::model_parameters(mods, ci =0.85) %>%select(Parameter, Coefficient, CI_low, CI_high) %>%filter(Parameter !="(Intercept)") %>%mutate(Model = mod_name,Informative =case_when( CI_low <0& CI_high <0~"yes", CI_low >0& CI_high >0~"yes",TRUE~"no")) return(ci_df) }
Consistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models (Sutherland et al. 2023). Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference (Arnold 2010).
Code
# get best models for all data TC_best_mods <-get.models(TC_selec, subset = delta <2, REML = T)# get best models for felidae subset dataTC_best_mods_F <-get.models(TC_selec_F, subset = delta <2, REML = T)# get best models for Canidae subset dataTC_best_mods_C <-get.models(TC_selec_C, subset = delta <2, REML = T)# get best models for Mustelidae subset dataTC_best_mods_M <-get.models(TC_selec_M, subset = delta <2)# Apply the function to obtain the table of coefficients of selected the modelsTC_best_ci <-map2_df(names(TC_best_mods), TC_best_mods, get_ci) TC_best_ci_F <-map2_df(names(TC_best_mods_F), TC_best_mods_F, get_ci)TC_best_ci_C <-map2_df(names(TC_best_mods_C), TC_best_mods_C, get_ci)TC_best_ci_M <-map2_df(names(TC_best_mods_M), TC_best_mods_M, get_cibeta) ci_plot(TC_best_ci)+labs(title="All data")ci_plot(TC_best_ci_F)+labs(title="Felidae dominant competitor")ci_plot(TC_best_ci_C)+labs(title="Canidae dominant competitor")ci_plot(TC_best_ci_M)+labs(title="Mustelidae dominant competitor")
All data best model 85% CI
Felidae data best model 85% CI
Canidae data best model 85% CI
Mustelidae data best model 85% CI
Summary selected models
Note that some variables have coefficients whose 85% confidence intervals overlap 0. These variables are retained when they have interactions with other variables and these interactions do not overlap with 0.
Arnold, Todd W. 2010. “Uninformative Parameters and Model Selection Using Akaike’s Information Criterion.”The Journal of Wildlife Management 74 (6): 1175–78. https://doi.org/10.2193/2009-367.
Brooks, Mollie E., Kasper Kristensen, Koen J. van, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Maechler, and Benjamin M. Bolker. 2017. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Modeling” 9. https://doi.org/10.32614/RJ-2017-066.
Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. New York: Springer-Verlag. //www.springer.com/us/book/9780387953649.
Faurby, Søren, Rasmus Ø Pedersen, Matt Davis, Simon D. Schowanek, Scott Jarvie, Alexandre Antonelli, and Jens-Christian Svenning. 2020. MegaPast2Future/PHYLACINE_1.2: PHYLACINE Version 1.2.1. Zenodo. https://doi.org/10.5281/zenodo.3690867.
Hassanin, Alexandre, Géraldine Veron, Anne Ropiquet, Bettine Jansen van Vuuren, Alexis Lécu, Steven M. Goodman, Jibran Haider, and Trung Thanh Nguyen. 2021. “Evolutionary History of Carnivora (Mammalia, Laurasiatheria) Inferred from Mitochondrial Genomes.”PLOS ONE 16 (2): e0240770. https://doi.org/10.1371/journal.pone.0240770.
Lüdecke, Daniel, Mattan S. Ben-Shachar, Indrajeet Patil, Philip Waggoner, and Dominique Makowski. 2021. “Performance: An r Package for Assessment, Comparison and Testing of Statistical Models” 6: 3139. https://doi.org/10.21105/joss.03139.
Makowski, Dominique, Mattan S. Ben-Shachar, Indrajeet Patil, and Daniel Lüdecke. 2020. “Estimation of Model-Based Predictions, Contrasts and Means.”https://github.com/easystats/modelbased.
Richmond, Orien M. W., James E. Hines, and Steven R. Beissinger. 2010. “Two-Species Occupancy Models: A New Parameterization Applied to Co-Occurrence of Secretive Rails.”Ecological Applications 20 (7): 2036–46. https://doi.org/10.1890/09-0470.1.
Ridout, M. S., and M. Linkie. 2009. “Estimating Overlap of Daily Activity Patterns from Camera Trap Data.”Journal of Agricultural, Biological, and Environmental Statistics 14 (3): 322–37. https://doi.org/10.1198/jabes.2009.08038.
Schloerke, Barret, Di Cook, Joseph Larmarange, Francois Briatte, Moritz Marbach, Edwin Thoen, Amos Elberg, and Jason Crowley. 2021. “GGally: Extension to ’Ggplot2’.”https://CRAN.R-project.org/package=GGally.
Soria, Carmen D., Michela Pacifici, Moreno Di Marco, Sarah M. Stephen, and Carlo Rondinini. 2021. “COMBINE: A Coalesced Mammal Database of Intrinsic and Extrinsic Traits.”Ecology 102 (6): e03344. https://doi.org/10.1002/ecy.3344.
Waddle, J.Hardin, Robert M. Dorazio, Susan C. Walls, Kenneth G. Rice, Jeff Beauchamp, Melinda J. Shuman, and Frank Mazzotti. 2010. “A New Parameterization for Estimating Co-Occurrence of Interacting Species.”Ecological Applications 20 (5): 302315. https://doi.org/https://doi.org/10.1890/09-0850.1.
Zuur, Alain F., Elena N. Ieno, and Chris S. Elphick. 2010. “A Protocol for Data Exploration to Avoid Common Statistical Problems.”Methods in Ecology and Evolution 1 (1): 3–14. https://doi.org/https://doi.org/10.1111/j.2041-210X.2009.00001.x.
Zuur, Alain F, Elena N Ieno, Neil J Walker, Anatoly A Saveliev, and Graham M Smith. 2009a. Mixed Effects Models and Extensions in Ecology with r. Vol. 574. Springer.
———. 2009b. Mixed Effects Models and Extensions in Ecology with r. 1st ed. Statistics for Biology and Health. Springer.
Source Code
---title: "Competitive interactions co-occurrence modelling"format: html: message: false warning: false code-fold: true code-tools: true toc: true self-contained: true theme: journalbibliography: references.bib---# Packages```{r}#| message: false#| warning: false# Data manipulation packagelibrary(tidyverse) # Easily Install and Load the 'Tidyverse'library(psych) # Procedures for Psychological, Psychometric, and Personality Research# Modelling packageslibrary(nlme) # Linear and Nonlinear Mixed Effects Modelslibrary(betareg) # Beta Regressionlibrary(glmmTMB) # Generalized Linear Mixed Models using Template Model Builderlibrary(AICcmodavg) # Model Selection and Multimodel Inference Based on (Q)AIC(c)library(MuMIn) # Multi-Model Inference#Modelling utility packageslibrary(broom) # Convert Statistical Objects into Tidy Tibbleslibrary(broom.mixed) # Tidying Methods for Mixed Modelslibrary(insight) # Easy Access to Model Information for Various Model Objectslibrary(modelbased) # Estimation of Model-Based Predictions, Contrasts and Meanslibrary(parameters) # Processing of Model Parameterslibrary(performance) # Assessment of Regression Models Performancelibrary(DHARMa) # Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Modelslibrary(ggeffects) # Create Tidy Data Frames of Marginal Effects for 'ggplot' from Model Outputs# Grapich packageslibrary(GGally) # Extension to 'ggplot2'library(patchwork) # The Composer of Plots# Table output packageslibrary(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntaxlibrary(DT)library(showtext)font_add_google("Roboto", "Roboto")showtext_auto()showtext_opts(dpi =300)```# Variables description| Variable | Description | Units | Source ||------------------|--------------------|------------------|------------------|| SIF | Spatial species interaction factor. SIF\<1: segregation, SIF=1: independence, SIF\>1= aggregation | index | Each study || Ov_coeff | Temporal overlap coefficient. Ov_coeff =1: complete activity overlap, Ov_Coeff =0: No activity overlap | index | Each study || CBMratio | Competitors body mass ratio: ln( dominant competitor body mass/ subordinate competitor body mass) | index | COMBINE data base [@soria2021] || Diet_dist | Dominant and subordinate competitor category diets: Same if pair of competitors shares de diet category, different if not | category | PHYLACINE trait data base [@faurby2020] || Phy_dist | Pairwise mitochondrial DNA phylogenetic distances | mitochondrial DNA phylogenetic distances | @hassanin2021 || Abs_lat | Absolute latitude of each study when reported | coordinates | Each study || Avg_dist | Average distance (km) between camera trap stations reported. Only used for spatial co-occurrence models | Average liner distance in km | Each study || Samp_dur | duration of sampling in months. Only used for temporal co-occurrence models | Months | Each study || Locality | Each sampling location where species co-occurrence information was obtained, within each study | category | Each study || Label | Each study unique ID | category | Each study |: List of variables used to modelling spatial co-occurrence patterns of competitive interactions {#tbl-variables}# DataThe data correspond to studies evaluating the spatial and temporal co-occurrence of mammals of the order Carnivora. In all cases the species are of the same order and present ecological interactions of competition.::: panel-tabset## Spatial dataSpatial co-occurrence data of carnivorous mammals were measured using the species interaction factor (SIF). This is a parameter derived from multi-species occupancy models [@waddle2010; @richmond2010].```{r}# Load the spatial co-ocurrence data base of competitive interactionsspatC_db <-read_delim("Data/Model data/spatC_db.csv", delim =";" )[,-1] %>%# Select the variablesselect(SIF, D_competitor, D_family, S_competitor, S_family, mass_ratio, diet_dist, Lat, Lat_abs, p_distance, Avg_dist, Locality, Label)DT::datatable(spatC_db)```The @tbl-Ssum contain the general description of the numeric variables of spatial data base. The table was constructed with psych package [@psych]```{r}#| label: tbl-Ssum#| tbl-cap: General description of the spatial co-occurrence data base for competitive interactions. Spatnum_summary <- spatC_db %>%select_if(is.numeric) %>%describe(. ,fast = T) kbl(Spatnum_summary, caption ="General description of the spatial co-occurrence data base for competitive interactions", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```We use the GGAlly package to explore the visual relationship between numerical variables [@GGally]. We used pearson's correlation coefficient```{r}#| message: false#| warning: false#| label: fig-Sccorplot#| fig-cap: Spatial competitive co-occurrence corplot(Spat_cor <-select_if(spatC_db, is.numeric) %>%# select numeric variablesggpairs(.,# Correlation coefficient upper partupper =list(continuous=wrap("cor", method="pearson", digits=2, corSize=80)),lower =list( continuous="smooth")) +theme_bw()+theme(text=element_text(family ="Roboto")))```## Temporal dataThe temporal data corresponds to the co-occurrence of Carnivora mammals measured by overlap coefficient of kernel activity curves [@ridout2009].```{r}tempC_db <-read_delim("Data/Model data/tempC_db.csv", delim =";")[,-1] %>%select(Ov_coeff, D_competitor, D_family, S_competitor, S_family, mass_ratio, diet_dist, Lat, Lat_abs, p_distance, Samp_dur, Locality, Label)DT::datatable(tempC_db )```The @tbl-Tsum contain the general description of the numeric variables of spatial data base.```{r}#| label: tbl-Tsum#| tbl-cap: General description of the temporal co-occurrence data base for competitive interactions.Tempnum_summary <- tempC_db %>%select_if(is.numeric) %>%describe(. ,fast = T) kbl(Tempnum_summary, caption ="General description of the temporal co-occurrence data base for competitive interactions", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```We use the GGAlly package to explore the visual relationship between numerical variables [@GGally]. We used pearson's correlation coefficient```{r}#| message: false#| warning: false#| label: fig-Scorplot#| fig-cap: Temporal competitive co-occurrence corplot(Temp_cor <-select_if(tempC_db, is.numeric) %>%ggpairs(.,upper =list(continuous=wrap("cor", method="pearson", digits=2, corSize=80)),lower =list( continuous="smooth")) +theme_bw()+theme(text=element_text(family ="Roboto")))```:::# Modelling procedure## Standardize numeric covariatesThe numerical covariates were standardized to mean 0 and standard deviation 1, which facilitates model fitting and direct comparison of regression coefficients. In addition, pairs of species that do not contain information on any of the variables (cells with NAs) were eliminated.Because we aimed to evaluate whether the patterns found varied depending on the dominant competitor family, we made subsets of the databases. For spatial co-occurrence we only subset the felidae family. For temporal co-occurrence we subset the families Felidae, Canidae and Mustelidae.```{r}# Function to standardized the numeric variablesscale_this <-function(data){scale(data) %>%as.vector() }# Select the variables for Spatial and temporal data setsSC_vars <-c("mass_ratio", "Lat_abs", "p_distance", "Avg_dist")TC_vars <-c("mass_ratio", "Lat_abs", "p_distance", "Samp_dur")# All data SpatialspatC_modsdf <- spatC_db %>%drop_na() %>%mutate(across(all_of(SC_vars), scale_this)) dim(spatC_modsdf)#Subset for spatial data# Felidae SpatialspatC_modsdf_F <- spatC_modsdf %>%filter(D_family =="Felidae")dim(spatC_modsdf_F)# All data TemporaltempC_modsdf <- tempC_db %>%drop_na() %>%mutate(across(all_of(TC_vars), scale_this))dim(tempC_modsdf)#subset for temporal data# Felidae TemporaltempC_modsdf_F <- tempC_modsdf %>%filter(D_family =="Felidae")dim(tempC_modsdf_F)# Canidae TemporaltempC_modsdf_C <- tempC_modsdf %>%filter(D_family =="Canidae")dim(tempC_modsdf_C)# Mustelidae TemporaltempC_modsdf_M <- tempC_modsdf %>%filter(D_family =="Mustelidae")dim(tempC_modsdf_M)```## Spatial modelling procedure### Influential observationsIdentify the presence of outliers or observations that may influence the models. For the spatial co-occurrence data we fit a general linear model with all covariates and interactions. We then checked for the presence of extreme data using Laverage distance plots from performance package [@performance].```{r}#| message: false#| warning: false#| fig-cap: Spatial co-occurrence leverage distance#| fig-subcap: #| - "Outliers of all data spatial model"#| - "Outliers of felidae dominant competior spatial model"#| layout-ncol: 2#| fig-width: 8# Fit model with all dataSCmods_outl <-lm(SIF ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2, data= spatC_modsdf)# Fit model with only felids as dominantSCmods_outl_F <-lm(SIF ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2, data= spatC_modsdf_F)# Check outliers for each modelSC_out <-check_outliers(SCmods_outl)SC_out_F <-check_outliers(SCmods_outl_F)# Leverage plot for all data modelplot(SC_out)+theme_bw()+theme(text=element_text(family ="Roboto"))# Leverage plot for felidae as dominant competitorplot(SC_out_F)+theme_bw()+theme(text=element_text(family ="Roboto"))```### CollinearityWe checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value \> 10 indicates a high collinearity of the variables [@zuur2010a].To calculate the VIF we use the Performance package [@performance].```{r}#| message: false#| warning: false#| label: fig-Scoll#| fig-cap: Spatial co-occurrence collinearity#| fig-subcap: #| - "VIF of all data spatial model"#| - "VIF of felidae dominant competior spatial model"#| layout-ncol: 2#| fig-width: 8# Fit model with all dataspatC_coll <-lm(SIF ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist, data= spatC_modsdf)# Fit model with only felids as dominantspatC_coll_F <-lm(SIF ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist, data= spatC_modsdf_F)# Check outliers for each modelspatC_colltab <-check_collinearity(spatC_coll)spatC_colltab_F <-check_collinearity(spatC_coll_F)# Vif plotsplot(spatC_colltab)+labs(subtitle ="All data Spatial model")+theme_bw()+theme(legend.position ="none",text=element_text(family ="Roboto"))plot(spatC_colltab_F)+labs(subtitle ="Felidae dominant competitor Spatial model")+theme_bw()+theme(legend.position ="none",text=element_text(family ="Roboto"))```### Model assumptions#### All dataWe first check the assumptions of the model. To do this we fit the most complex model containing all available variable interactions and random factors. We then visually evaluated the normality of the residuals, linearity and homogeneity of variance using the performance package [@performance].```{r}Spat_diag_m <-lme(SIF ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs)^2, random =~1|Label/Locality,method ="REML",data= spatC_modsdf)Spat_diagnostic <-check_model(Spat_diag_m, check =c("homogeneity", "linearity", "qq", "normality", "reqq"))Spat_diagnostic ```We detect deviations from assumptions. Particularly there is a kurtosis in the distribution of the residuals. There are also deviations from linearity. For that reason, we fit a T-student family glm using glmmTMB package [@glmmTMB]. To evaluate the GLMM we use the DHARMa package [@DHARMa], which applies simulated residuals to evaluate the assumptions.```{r}Spat_diag_t <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label/Locality),data=spatC_modsdf, family=t_family(link ="identity"), REML = T)res_t <-simulateResiduals(Spat_diag_t)plot(res_t)```The t-family GLMM improves considerably and the graphs suggest that it meets the distribution assumptions.#### FelidaeThis section verifies the model assumptions considering only felids as dominant competitors.```{r}Spat_diag_m_F <-lme(SIF ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2, random =~1|Label/Locality,method ="REML",data= spatC_modsdf_F)Spat_diagnostic_F <-check_model(Spat_diag_m_F, check =c("homogeneity", "linearity", "qq", "normality"))Spat_diagnostic_F```Deviations from the assumptions are detected, so we tried a GLMM of family t student.```{r}Spat_diag_tF <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label/Locality),data=spatC_modsdf_F, family=t_family(link ="identity"), REML = T)res_t_F <-simulateResiduals(Spat_diag_tF)plot(res_t_F)```Since deviations from homogeneity of variance were detected, we evaluated which variables are involved in the model inadequacies.```{r}#| message: false#| warning: false#| fig-subcap: #| - "Uniformity of variance test of Diet category"#| - "Uniformity of variance test of mass ratio"#| - "Uniformity of variance test of absolute latitude"#| - "Uniformity of variance test of phylogenetic distance"#| - "Uniformity of variance test of Average distance"#| layout-ncol: 2#| fig-width: 8testCategorical(res_t_F, spatC_modsdf_F$diet_dist) testQuantiles(res_t_F, spatC_modsdf_F$mass_ratio)testQuantiles(res_t_F, spatC_modsdf_F$Lat_abs)testQuantiles(res_t_F, spatC_modsdf_F$p_distance)testQuantiles(res_t_F, spatC_modsdf_F$Avg_dist)```We detected non-uniformity in the range of body mass ratio, phylogenetic distance, absolute latitude and average camera distance. Additionally, we identified that observation 84 (previously identified as outlier) does have an effect on model fit. To improve the fit, we modeled the variation of the variables by means of the dispformula term, from the glmm TMB package [@glmmTMB].```{r}Spat_diag_tF2 <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label/Locality),dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F[-84,], family=t_family(link ="identity"), REML = T)res_t_F2 <-simulateResiduals(Spat_diag_tF2)plot(res_t_F2)```### Random structure#### All dataFollowing the @zuur2009 protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the previously selected model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package [@AICcmodavg].```{r}#| tbl-cap: Random-effects structure selection table for spatial co-occurrence data# Fit model without random effectsSpat_r0 <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2,data=spatC_modsdf, family=t_family(link ="identity"), REML = T)# Fit model with label as random effectSpat_r1 <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label),data=spatC_modsdf, family=t_family(link ="identity"), REML = T)# Fit model with label and Locality as random effectsSpat_r2 <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label/Locality),data=spatC_modsdf, family=t_family(link ="identity"), REML = T)# Rank the models with the AICspat_AICtab <-aictab(cand.set =list(Spat_r0, Spat_r1, Spat_r2),modnames =c("no random","Label random", "Label/Locality random"),second.ord = F)kbl(spat_AICtab, caption ="Random-effects structure selection table for spatial co-occurrence data", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```The @tbl-Sraic suggests that both models with Label as random effects and without random effects are equally supported. To visualize the variation of the random parameters of each group, we use the estimate_grouplevel function of the modelbased package [@modelbased].```{r}estimate_grouplevel(Spat_r1) %>%# Get random parametersplot()+#plot# aesthetic changestheme_bw()+theme(text=element_text(family ="Roboto"))```There is no significant variation from each group level in random effects. We now check the goodness of fit of these models.```{r}#| message: false#| warning: false#| fig-subcap: #| - "Gof-test of model with Label as random effects"#| - "Gof-test of model without random effects"#| layout-ncol: 2#| fig-width: 8simulateResiduals(Spat_r1, plot = T)simulateResiduals(Spat_r0, plot = T)```Both models have no lack of fit. In this case, we select the simpler one (without random effects) to continue with the selection of fixed effects.#### FelidaeWe performed the same procedure for Felidae```{r}Spat_r0_F <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2,dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F[-84,], family=t_family(link ="identity"), REML = T,start =list(psi =log(2.87)),map =list(psi =factor(NA)))Spat_r1_F <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label),dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F[-84,], family=t_family(link ="identity"), REML = T,start =list(psi =log(2.87)),map =list(psi =factor(NA)))Spat_r2_F <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2+(1|Label/Locality),dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F[-84,], family=t_family(link ="identity"), REML = T,start =list(psi =log(2.87)),map =list(psi =factor(NA)))spat_AICtab_F <-aictab(cand.set =list(Spat_r0_F, Spat_r1_F, Spat_r2_F),modnames =c("no random","Label random", "Label/Locality random"),second.ord = F)kbl(spat_AICtab_F, caption ="Random-effects structure selection table for spatial co-occurrence data when Felidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```In the case of Felidae the model selection indicates that random effects are not adequate. This is because there is no variation between Label or Locality levels.```{r}estimate_grouplevel(Spat_r1_F) %>%plot()+theme_bw()+theme(text=element_text(family ="Roboto"))```All label coefficient are the same for felids data### Fixed predictor selectionsTo select the fixed effects structure we also used the Akaike information criterion (AIC). Models with a delta AIC \<2 are considered equally plausible [@burnham2002]. Because we are going to choose the fixed effects structure, we refit the previously identified model without using the restricted maximum likelihood (REML) [@zuur2009a]. Since co-occurrence patterns may be generated by the interaction of variables, we fit all possible combinations of variables, limiting to a maximum of three parameters per model. To do this we will used the dredge function of the MuMIn package [@MuMIn].#### All data```{r}# Create a global modelSC_global <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2,data=spatC_modsdf, family=t_family(link ="identity"), REML = F,na.action ="na.fail",start =list(psi =log(1.69)),map =list(psi =factor(NA)))#Fit all posible combination limit to 3 parametersSC_selec <-dredge(SC_global, rank ="AIC", m.lim=c(NA,3))``````{r}kbl(SC_selec, caption ="Fixed-effects structure selection table for spatial co-occurrence data", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```#### FelidaeNote that this time the limit number of parameters in the function is 7. This is due to the number of parameters of the dispersion terms, which add 4 additional parameters.```{r}spatC_modsdf_F <- spatC_modsdf_F[-84,]SC_global_F <-glmmTMB(formula = SIF~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Avg_dist)^2,dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F, family=t_family(link ="identity"), REML = F,start =list(psi =log(2.87)),map =list(psi =factor(NA)),na.action ="na.fail")SC_selec_F <-dredge(SC_global_F, rank ="AIC", fixed =c("disp(Avg_dist)", "disp(Lat_abs)", "disp(mass_ratio)", "disp(p_distance)"), # Mantain dispersion parameters in all modelsm.lim=c(NA,7)) # 7 parameteres becasue we must add the dispersion terms``````{r}kbl(SC_selec_F, caption ="Fixed-effects structure selection table for spatial co-occurrence data when Felidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```### Confidence intervals explorationsConsistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models [@sutherland2023]. Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference [@arnold2010].```{r}# Function to get 85%ICget_ci <-function(mod_name, mods){ ci_df <- parameters::model_parameters(mods, ci =0.85) %>%select(Parameter, Coefficient, CI_low, CI_high, Component) %>%filter(Component =="conditional") %>%filter(Parameter !="(Intercept)") %>%mutate(Model = mod_name,Informative =case_when( CI_low <0& CI_high <0~"yes", CI_low >0& CI_high >0~"yes",TRUE~"no")) return(ci_df) }# Function to plotci_plot <-function(ci_df){ggplot(ci_df, aes(x=Coefficient, y= Parameter))+geom_pointrange(aes(xmin=CI_low, xmax= CI_high,col= Model, linetype= Informative),position =position_dodge2(0.5), linewidth=1)+scale_linetype_manual(values =c("no"="dashed", "yes"="solid"))+geom_vline(xintercept =0, linetype="dashed")+scale_color_viridis_d()+labs(caption ="estimates with 85% ci intervals",col="Model ID")+theme_bw()+theme(text=element_text(family ="Roboto"))}``````{r}#| message: false#| warning: false#| fig-subcap: #| - "All data best model 85% CI"#| - "Felidae data best model 85% CI"#| layout-ncol: 2#| fig-width: 8# get best models for all data SC_best_mods <-get.models(SC_selec, subset = delta <2)# get best models for felidae subset dataSC_best_mods_F <-get.models(SC_selec_F, subset = delta <2)# Apply the function to obtain the table of coefficients of selected the modelsSC_best_ci <-map2_df(names(SC_best_mods), SC_best_mods, get_ci) SC_best_ci_F <-map2_df(names(SC_best_mods_F), SC_best_mods_F, get_ci) ci_plot(SC_best_ci)+labs(title="All data")ci_plot(SC_best_ci_F)+labs(title="Felidae dominant competitor")```## Temporal modelling procedure### Influential observationsIdentify the presence of outliers or observations that may influence the models. For the temporal co-ocurrence data we fit a beta family generalized linear model with all covariates and interactions [@betareg]. We then checked for the presence of extreme observations using Cook's distance plots from performance package [@performance].```{r}#| message: false#| warning: false#| fig-cap: Temporal co-occurrence cook distance distance#| fig-subcap: #| - "Outliers of all data temporal model"#| - "Outliers of Felidae dominant competior Temporal model"#| - "Outliers of Celidae dominant competior Temporal model"#| - "Outliers of Mustelidae dominant competior Temporal model"#| layout-ncol: 2#| fig-width: 8# Models for each dataTCmods_outl <-betareg(Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2, data= tempC_modsdf)TCmods_outl_F <-betareg(Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2, data= tempC_modsdf_F)TCmods_outl_C <-betareg(Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2, data= tempC_modsdf_C)TCmods_outl_M <-betareg(Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2, data= tempC_modsdf_M)#Check outliersTC_out <-check_outliers(TCmods_outl)TC_out_F <-check_outliers(TCmods_outl_F)TC_out_C <-check_outliers(TCmods_outl_C)TC_out_M <-check_outliers(TCmods_outl_M)# Cook's distance plotsplot(TC_out)+theme_bw()+theme(text=element_text(family ="Roboto"))plot(TC_out_F)+theme_bw()+theme(text=element_text(family ="Roboto"))plot(TC_out_C)+theme_bw()+theme(text=element_text(family ="Roboto"))plot(TC_out_M)+theme_bw()+theme(text=element_text(family ="Roboto"))```A possible outlier was detected for the Mustelidae data as a dominant competitor. It corresponds to observation 27 which exceeded the 0.98 Cook's distance threshold.### CollinearityWe checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value \> 10 indicates a high collinearity of the variables [@zuur2010a].To calculate the VIF we use the Performance package [@performance].```{r}#| label: fig-TCcoll_fam#| fig-cap: Outliers model detection temporal models#| fig-subcap: #| - "VIF of all data model"#| - "VIF of Felidae model"#| - "VIF of Canidae model"#| - "VIF of Mustelidae model"#| layout-ncol: 2#| fig-width: 8# Fit modelstempC_coll <-betareg(Ov_coeff ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur, data= tempC_modsdf)tempC_coll_F <-betareg(Ov_coeff ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur, data= tempC_modsdf_F)tempC_coll_C <-betareg(Ov_coeff ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur, data= tempC_modsdf_C)tempC_coll_M <-betareg(Ov_coeff ~ mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur, data= tempC_modsdf_M)# Check collinearitytempC_colltab <-check_collinearity(tempC_coll)tempC_colltab_F <-check_collinearity(tempC_coll_F)tempC_colltab_C <-check_collinearity(tempC_coll_C)tempC_colltab_M <-check_collinearity(tempC_coll_M )# Vif plotsplot(tempC_colltab)+labs(subtitle ="All data Temporal model")+theme_bw()+theme(legend.position ="none",text=element_text(family ="Roboto"))plot(tempC_colltab_F)+labs(subtitle ="Felidae Temporal model")+theme_bw()+theme(legend.position ="none",text=element_text(family ="Roboto"))plot(tempC_colltab_C)+labs(subtitle ="Canidae Temporal model")+theme_bw()+theme(legend.position ="none",text=element_text(family ="Roboto"))plot(tempC_colltab_M)+labs(subtitle ="Mustelidae Temporal model")+theme_bw()+theme(legend.position ="none",text=element_text(family ="Roboto"))```### Model assumptionsWe verify the model assumptions by visual inspection of the simulated residuals from the DHARMa package [@DHARMa]. To do so we fit the more complex model which includes fixed variables and their interactions, as well as random effects.#### All data```{r}Temp_diag <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),data=tempC_modsdf, family=beta_family(), REML = T)Temp_diag_res <-simulateResiduals(Temp_diag)plot(Temp_diag_res)```Since deviations from homogeneity of variance were detected, we evaluated which variables are involved in the model inadequacies.```{r}#| message: false#| warning: false#| fig-subcap: #| - "Uniformity of variance test of Diet category"#| - "Uniformity of variance test of mass ratio"#| - "Uniformity of variance test of absolute latitude"#| - "Uniformity of variance test of phylogenetic distance"#| - "Uniformity of variance test of sampling duration"#| layout-ncol: 2#| fig-width: 8testCategorical(Temp_diag_res, tempC_modsdf$diet_dist) testQuantiles(Temp_diag_res , tempC_modsdf$mass_ratio)testQuantiles(Temp_diag_res , tempC_modsdf$Lat_abs)testQuantiles(Temp_diag_res , tempC_modsdf$p_distance)testQuantiles(Temp_diag_res , tempC_modsdf$Samp_dur)```We detected non-uniformity in the range of diet distance, phylogenetic distance, absolute latitude and mean camera distance. Additionally, we identified that observation 84 (previously identified as outlier) does have an effect on model fit. To improve the fit, we modeled the variation of the variables by means of the dispformula term, from the glmm TMB package [@glmmTMB].```{r}Temp_diag2 <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),dispformula =~diet_dist+ p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = T)Temp_diag_res2 <-simulateResiduals(Temp_diag2)plot(Temp_diag_res2)```#### Felidae```{r}Temp_diag_F <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),data=tempC_modsdf_F, family=beta_family(), REML = T)Temp_diag_res_F <-simulateResiduals(Temp_diag_F)plot(Temp_diag_res_F)``````{r}#| message: false#| warning: false#| fig-subcap: #| - "Uniformity of variance test of Diet category"#| - "Uniformity of variance test of mass ratio"#| - "Uniformity of variance test of absolute latitude"#| - "Uniformity of variance test of phylogenetic distance"#| - "Uniformity of variance test of Sampling duration"#| layout-ncol: 2#| fig-width: 8testCategorical(Temp_diag_F, tempC_modsdf_F$diet_dist)testQuantiles(Temp_diag_F, tempC_modsdf_F$mass_ratio)testQuantiles(Temp_diag_F, tempC_modsdf_F$Lat_abs)testQuantiles(Temp_diag_F, tempC_modsdf_F$p_distance)testQuantiles(Temp_diag_F, tempC_modsdf_F$Samp_dur)``````{r}Temp_diag_F2 <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),dispformula =~mass_ratio+p_distance+Lat_abs+Samp_dur,data=tempC_modsdf_F, family=beta_family(), REML = T)Temp_diag_res_F2 <-simulateResiduals(Temp_diag_F2)plot(Temp_diag_res_F2)```#### Canidae```{r}Temp_diag_C <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),data=tempC_modsdf_C, family=beta_family(), REML = T)Temp_diag_res_C <-simulateResiduals(Temp_diag_C)plot(Temp_diag_res_C)``````{r}#| message: false#| warning: false#| fig-subcap: #| - "Uniformity of variance test of Diet category"#| - "Uniformity of variance test of mass ratio"#| - "Uniformity of variance test of absolute latitude"#| - "Uniformity of variance test of phylogenetic distance"#| - "Uniformity of variance test of Sampling duration"#| layout-ncol: 2#| fig-width: 8testCategorical(Temp_diag_C, tempC_modsdf_C$diet_dist)testQuantiles(Temp_diag_C, tempC_modsdf_C$mass_ratio)testQuantiles(Temp_diag_C, tempC_modsdf_C$Lat_abs)testQuantiles(Temp_diag_C, tempC_modsdf_C$p_distance)testQuantiles(Temp_diag_C, tempC_modsdf_C$Samp_dur)``````{r}Temp_diag_C2 <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),dispformula =~Samp_dur+Lat_abs,data=tempC_modsdf_C, family=beta_family(), REML = T)Temp_diag_res_C2 <-simulateResiduals(Temp_diag_C2)plot(Temp_diag_res_C2)```#### Mustelidae```{r}Temp_diag_M <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),data=tempC_modsdf_M, family=beta_family(), REML = T)Temp_diag_res_M <-simulateResiduals(Temp_diag_M)plot(Temp_diag_res_M)```The previously identified outlier caused the inadequacies of the model.```{r}Temp_diag_M2 <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),data=tempC_modsdf_M[-27,], family=beta_family(), REML = T)Temp_diag_res_M2 <-simulateResiduals(Temp_diag_M2)plot(Temp_diag_res_M2)```### Random structureFollowing the @zuur2009 protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the previously selected model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package [@AICcmodavg].#### All data```{r}Temp_r0 <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2,dispformula =~diet_dist+ p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = T)Temp_r1 <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label),dispformula =~diet_dist+ p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = T)Temp_r2 <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label/Locality),dispformula =~diet_dist+ p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = T)temp_AICtab <-aictab(cand.set =list(Temp_r0, Temp_r1, Temp_r2),modnames =c("no random","Label random", "Label/Locality random"),second.ord = F)kbl(temp_AICtab, caption ="Random-effects structure selection table for Temporal co-occurrence data", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)``````{r}estimate_grouplevel(Temp_r1) %>%# Get random parametersplot()+#plot# aesthetic changestheme_bw()+theme(text=element_text(family ="Roboto"))```#### Felidae```{r}Temp_r0_F <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2,data=tempC_modsdf_F, family=beta_family(), REML = T)Temp_r1_F <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label),data=tempC_modsdf_F, family=beta_family(), REML = T)Temp_r2_F <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label/Locality),data=tempC_modsdf_F, family=beta_family(), REML = T)temp_AICtab_F <-aictab(cand.set =list(Temp_r0_F, Temp_r1_F, Temp_r2_F),modnames =c("no random","Label random", "Label/Locality random"),second.ord = F)kbl(temp_AICtab_F, caption ="Random-effects structure selection table for Temporal co-occurrence data when Felidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)``````{r}estimate_grouplevel(Temp_r1_F) %>%# Get random parametersplot()+#plot# aesthetic changestheme_bw()+theme(text=element_text(family ="Roboto"))```#### Canidae```{r}Temp_r0_C <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2,data=tempC_modsdf_C, family=beta_family(), REML = T)Temp_r1_C <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label),data=tempC_modsdf_C, family=beta_family(), REML = T)Temp_r2_C <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label/Locality),data=tempC_modsdf_C, family=beta_family(), REML = T)temp_AICtab_C <-aictab(cand.set =list(Temp_r0_C, Temp_r1_C, Temp_r2_C),modnames =c("no random","Label random", "Label/Locality random"),second.ord = F)kbl(temp_AICtab_C, caption ="Random-effects structure selection table for Temporal co-occurrence data when Canidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)``````{r}estimate_grouplevel(Temp_r1_C) %>%# Get random parametersplot()+#plot# aesthetic changestheme_bw()+theme(text=element_text(family ="Roboto"))```#### Mustelidae```{r}Temp_r0_M <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2,data=tempC_modsdf_M, family=beta_family(), REML = T)Temp_r1_M <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label),data=tempC_modsdf_M, family=beta_family(), REML = T)Temp_r2_M <-glmmTMB(formula = Ov_coeff~ (mass_ratio+I(mass_ratio^2)+ p_distance+ diet_dist+ Lat_abs+ Samp_dur)^2+(1|Label/Locality),data=tempC_modsdf_M, family=beta_family(), REML = T)temp_AICtab_M <-aictab(cand.set =list(Temp_r0_M, Temp_r1_M, Temp_r2_M),modnames =c("no random","Label random", "Label/Locality random"),second.ord = F)kbl(temp_AICtab_M, caption ="Random-effects structure selection table for Temporal co-occurrence data when Mustelidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)``````{r}estimate_grouplevel(Temp_r1_M) %>%# Get random parametersplot()+#plot# aesthetic changestheme_bw()+theme(text=element_text(family ="Roboto"))```### Fixed predictor selections#### All data```{r}TC_global <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label),dispformula =~p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = F,na.action ="na.fail")TC_selec <-dredge(TC_global, rank ="AIC", fixed =c("disp(p_distance)","disp(Lat_abs)","disp(Samp_dur)"), m.lim=c(NA,6))``````{r}kbl(TC_selec, caption ="Fixed-effects structure selection table for temporal co-occurrence data", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```#### Felidae```{r}TC_global_F <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label),dispformula =~p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf_F, family=beta_family(), REML = F,na.action ="na.fail")TC_selec_F <-dredge(TC_global_F, rank ="AIC", fixed =c("disp(p_distance)","disp(Lat_abs)","disp(Samp_dur)"), m.lim=c(NA,6))``````{r}kbl(TC_selec_F, caption ="Fixed-effects structure selection table for temporal co-occurrence data when Felidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```#### Canidae```{r}TC_global_C <-glmmTMB(formula = Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label),dispformula =~Samp_dur+Lat_abs,data=tempC_modsdf_C, family=beta_family(), REML = F,na.action ="na.fail")TC_selec_C <-dredge(TC_global_C, rank ="AIC", fixed =c("disp(Lat_abs)","disp(Samp_dur)"), m.lim=c(NA,5))``````{r}kbl(TC_selec_C, caption ="Fixed-effects structure selection table for temporal co-occurrence data when Canidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```#### Mustealidae```{r}tempC_modsdf_M2 <-slice(tempC_modsdf_M, -27)TC_global_M <-betareg(Ov_coeff ~ (mass_ratio+I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2,data=tempC_modsdf_M2, na.action ="na.fail")TC_selec_M <-dredge(TC_global_M, rank ="AIC", m.lim=c(NA,3))``````{r}kbl(TC_selec_M, caption ="Fixed-effects structure selection table for temporal co-occurrence data when Mustelidae is dominant competitor", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F)```### Confidence intervals explorations```{r}# Function to get 85%ICget_cibeta <-function(mod_name, mods){ ci_df <- parameters::model_parameters(mods, ci =0.85) %>%select(Parameter, Coefficient, CI_low, CI_high) %>%filter(Parameter !="(Intercept)") %>%mutate(Model = mod_name,Informative =case_when( CI_low <0& CI_high <0~"yes", CI_low >0& CI_high >0~"yes",TRUE~"no")) return(ci_df) }```Consistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models [@sutherland2023]. Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference [@arnold2010].```{r}#| message: false#| warning: false#| fig-subcap: #| - "All data best model 85% CI"#| - "Felidae data best model 85% CI"#| - "Canidae data best model 85% CI"#| - "Mustelidae data best model 85% CI"#| layout-ncol: 2#| fig-width: 8# get best models for all data TC_best_mods <-get.models(TC_selec, subset = delta <2, REML = T)# get best models for felidae subset dataTC_best_mods_F <-get.models(TC_selec_F, subset = delta <2, REML = T)# get best models for Canidae subset dataTC_best_mods_C <-get.models(TC_selec_C, subset = delta <2, REML = T)# get best models for Mustelidae subset dataTC_best_mods_M <-get.models(TC_selec_M, subset = delta <2)# Apply the function to obtain the table of coefficients of selected the modelsTC_best_ci <-map2_df(names(TC_best_mods), TC_best_mods, get_ci) TC_best_ci_F <-map2_df(names(TC_best_mods_F), TC_best_mods_F, get_ci)TC_best_ci_C <-map2_df(names(TC_best_mods_C), TC_best_mods_C, get_ci)TC_best_ci_M <-map2_df(names(TC_best_mods_M), TC_best_mods_M, get_cibeta) ci_plot(TC_best_ci)+labs(title="All data")ci_plot(TC_best_ci_F)+labs(title="Felidae dominant competitor")ci_plot(TC_best_ci_C)+labs(title="Canidae dominant competitor")ci_plot(TC_best_ci_M)+labs(title="Mustelidae dominant competitor")```## Summary selected modelsNote that some variables have coefficients whose 85% confidence intervals overlap 0. These variables are retained when they have interactions with other variables and these interactions do not overlap with 0.### Spatial co-occurrence models```{r}SC_mod1 <-glmmTMB(formula = SIF~I(mass_ratio^2)+ p_distance,data=spatC_modsdf, family=t_family(link ="identity"),REML = F, start =list(psi =log(1.69)),map =list(psi =factor(NA))) SC_mod2 <-glmmTMB(formula = SIF~I(mass_ratio^2),data=spatC_modsdf, family=t_family(link ="identity"),REML = F, start =list(psi =log(1.69)), map =list(psi =factor(NA))) SC_mods <-list(SC_mod1, SC_mod2) %>%map(model_parameters, digits=2, ci=0.85) %>%reduce(rbind) %>%mutate(Model=c(rep("SIF~ I(mass_ratio^2)+ p_distance", 4),rep("SIF~ I(mass_ratio^2)", 3))) %>%select(-SE, -z, -df_error) ``````{r}SC_mod1_F <-glmmTMB(formula = SIF~I(mass_ratio^2)*Lat_abs,dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F,family=t_family(link ="identity"),REML = F,start =list(psi =log(2.87)),map =list(psi =factor(NA))) SC_mod2_F <-glmmTMB(formula = SIF~I(mass_ratio^2)+ Lat_abs,dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F,family=t_family(link ="identity"),REML = F,start =list(psi =log(2.87)),map =list(psi =factor(NA))) SC_mod3_F <-glmmTMB(formula = SIF~I(mass_ratio^2),dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F,family=t_family(link ="identity"),REML = F,start =list(psi =log(2.87)),map =list(psi =factor(NA)))SC_mod4_F <-glmmTMB(formula = SIF~I(mass_ratio^2)*p_distance,dispformula =~ mass_ratio+p_distance+Lat_abs+ Avg_dist,data=spatC_modsdf_F,family=t_family(link ="identity"), REML = F,start =list(psi =log(2.87)), map =list(psi =factor(NA))) SC_mods_F <-list(SC_mod1_F, SC_mod2_F, SC_mod3_F, SC_mod4_F) %>%map(model_parameters, digits=2, ci=0.85) %>%reduce(rbind) %>%mutate(Model=c(rep("SIF~ I(mass_ratio^2)* Lat_abs", 9),rep("SIF~ I(mass_ratio^2)+ Lat_abs", 8),rep("SIF~ I(mass_ratio^2)", 7),rep("SIF~ I(mass_ratio^2)*p_distance", 9))) %>%select(-SE, -z, -df_error)``````{r}SC_mods_tab <-rbind(SC_mods, SC_mods_F)kbl(SC_mods_tab, caption ="Spatial co-occurrence selected models", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F) %>%pack_rows("Spatial co-occurrence models -all data",1, nrow(SC_mods)) %>%pack_rows("Spatial co-occurrence models -Felidae as dominant competitor",1+nrow(SC_mods),nrow(SC_mods)+nrow(SC_mods_F))```#### Goftest of spatial selected models```{r}#| message: false#| warning: false#| fig-subcap: #| - "SIF~ I(mass_ratio^2)+ p_distance goftest -all data"#| - "SIF~ I(mass_ratio^2) goftest - -all data"#| - "SIF~ p_distance goftest- -all data"#| layout-ncol: 2#| fig-width: 8SC_best_goft_list <-list(SC_mod1, SC_mod2) SC_best_goft <-lapply(SC_best_goft_list, simulateResiduals, plot=T)``````{r}#| message: false#| warning: false#| fig-subcap: #| - "SSIF~ I(mass_ratio^2)* Lat_abs goftest -Felidae as dominant competitor"#| - "SIF~ I(mass_ratio^2)+ Lat_abs goftest -Felidae as dominant competitor"#| - "SIF~ I(mass_ratio^2) goftest -Felidae as dominant competitor"#| - "SIF~ I(mass_ratio^2)*p_distance goftest -Felidae as dominant competitor"#| layout-ncol: 2#| fig-width: 8SC_best_goft_list_F <-list(SC_mod1_F, SC_mod2_F, SC_mod3_F, SC_mod4_F) SC_best_goft_F <-lapply(SC_best_goft_list_F, simulateResiduals, plot=T)```Because all models fit adequately, we base our inference on models that contain interactions of the variables### Temporal co-occurrence models```{r}TC_mod1 <-glmmTMB(formula = Ov_coeff ~ Lat_abs+ mass_ratio+ p_distance + (1|Label),dispformula =~p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = T)TC_mod2 <-glmmTMB(formula = Ov_coeff ~ Lat_abs+I(mass_ratio^2)+ p_distance+ (1|Label),dispformula =~p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = T)TC_mod3 <-glmmTMB(formula = Ov_coeff ~I(mass_ratio^2)*p_distance+ (1|Label),dispformula =~p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf, family=beta_family(), REML = T)TC_mods <-list(TC_mod1, TC_mod2, TC_mod3) %>%map(model_parameters, digits=2, ci=0.85) %>%reduce(rbind) %>%mutate(Model=c(rep("Ov_coeff ~ Lat_abs+ mass_ratio+ p_distance + (1|Label)", 9),rep(" Ov_coeff ~ Lat_abs+ I(mass_ratio^2)+ p_distance+ (1|Label)", 9),rep("Ov_coeff ~ I(mass_ratio^2)*p_distance+ (1|Label)", 9))) %>%select(-SE, -z, -df_error, -Group)``````{r}TC_mod1_F <-glmmTMB(formula = Ov_coeff ~ p_distance*diet_dist+ (1|Label),dispformula =~p_distance+ Samp_dur+Lat_abs,data=tempC_modsdf_F, family=beta_family(), REML = T)TC_mods_F <-model_parameters(TC_mod1_F, digits=2, ci=0.85) %>%mutate(Model=c(rep("Ov_coeff ~ p_distance*diet_dist+ (1|Label)", 9))) %>%select(-SE, -z, -df_error, -Group)``````{r}TC_mod1_C <-glmmTMB(formula = Ov_coeff ~ Lat_abs+ mass_ratio+ (1|Label),dispformula =~Samp_dur+Lat_abs,data=tempC_modsdf_C, family=beta_family(), REML = T)TC_mods_C <-list(TC_mod1_C) %>%map(model_parameters, digits=2, ci=0.85) %>%reduce(rbind) %>%mutate(Model=c(rep("Ov_coeff ~ Lat_abs+ mass_ratio+ (1|Label)", 7))) %>%select(-SE, -z, -df_error,-Group)``````{r}TC_mod1_M<-betareg(Ov_coeff ~ p_distance + diet_dist+Lat_abs,data=tempC_modsdf_M2)TC_mods_M <-model_parameters(TC_mod1_M, digits=2, ci=0.85) %>%mutate(Component="conditional",Effects="fixed",Model="Ov_coeff ~ p_distance + diet_dist+Lat_abs")%>%select(-SE, -z, -df_error)``````{r}TC_mods_tab <-rbind(TC_mods, TC_mods_F, TC_mods_C, TC_mods_M)kbl(TC_mods_tab, caption ="Temporal co-occurrence selected models", digits =3) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width = F) %>%pack_rows("Temporal co-occurrence models -all data",1, nrow(TC_mods)) %>%pack_rows("Temporal co-occurrence models -Felidae as dominant competitor",1+nrow(TC_mods),nrow(TC_mods)+nrow(TC_mods_F)) %>%pack_rows("Temporal co-occurrence models -Canidae as dominant competitor",1+nrow(TC_mods)+nrow(TC_mods_F), nrow(TC_mods)+nrow(TC_mods_F)+nrow(TC_mods_C)) %>%pack_rows("Temporal co-occurrence models -Mustelidae as dominant competitor", 1+nrow(TC_mods)+nrow(TC_mods_F)+nrow(TC_mods_C), nrow(TC_mods)+nrow(TC_mods_F)+nrow(TC_mods_C)+nrow(TC_mods_M))```#### Goftest of temporal selected models```{r}#| message: false#| warning: false#| fig-subcap: #| - "Ov_coeff ~ Lat_abs+ mass_ratio+ p_distance + (1|Label) goftest -all data"#| - "Ov_coeff ~ Lat_abs+ I(mass_ratio^2)+ p_distance+ (1|Label) goftest - -all data"#| - "Ov_coeff ~ I(mass_ratio^2)*p_distance+ (1|Label) goftest- -all data"#| layout-ncol: 2#| fig-width: 8TC_best_goft_list <-list(TC_mod1, TC_mod2, TC_mod3) TC_best_goft <-lapply(TC_best_goft_list, simulateResiduals, plot=T)``````{r}#| message: false#| warning: false#| fig-subcap: #| - "Ov_coeff ~ p_distance*diet_dist+ (1|Label) goftest -Felidae as dominant competitor"#| layout-ncol: 2#| fig-width: 8TC_best_goft_list_F <-list(TC_mod1_F) TC_best_goft_F <-lapply(TC_best_goft_list_F, simulateResiduals, plot=T)``````{r}#| message: false#| warning: false#| fig-subcap: #| - "Ov_coeff ~ mass_ratio + (1 | Label) goftest -Canidae as dominant competitor"#| - "Ov_coeff ~ Lat_abs + mass_ratio + (1 | Label) goftest -Canidae as dominant competitor"#| layout-ncol: 2#| fig-width: 8TC_best_goft_list_C <-list(TC_mod1_C) TC_best_goft_C <-lapply(TC_best_goft_list_C, simulateResiduals, plot=T)``````{r}#| message: false#| warning: false#| layout-ncol: 2#| fig-width: 8plot(TC_mod1_M)```## Predictions```{r}# Spatial all dataSC_pred_mass <-ggemmeans(SC_mod1, terms =c("mass_ratio[-1.2:3.1 by=.2]"), ci_level=.85,ci.lvl = .85) %>%mutate(mass_real= (x*attr(scale(spatC_db$mass_ratio), "scaled:scale"))+attr(scale(spatC_db$mass_ratio), "scaled:center"))SC_pred_pdist <-ggemmeans(SC_mod1, terms ="p_distance[-2.3:1.2 by=.2]", ci_level=.85,ci.lvl = .85) %>%mutate(p_real= (x*attr(scale(spatC_db$p_distance), "scaled:scale"))+attr(scale(spatC_db$p_distance), "scaled:center"))# Spatial FelidaeSCp_distvec_F <-c(max(spatC_modsdf_F$p_distance), mean(spatC_modsdf_F$p_distance),min(spatC_modsdf_F$p_distance))SCp_latvec_F <-c(max(spatC_modsdf_F$Lat_abs), mean(spatC_modsdf_F$Lat_abs),min(spatC_modsdf_F$Lat_abs))SC_pred_massxp_F <-ggemmeans(SC_mod4_F, terms =c("mass_ratio[all]", "p_distance[SCp_distvec_F]"),ci_level=.85,ci.lvl = .85) %>%mutate( mass_real= (x*attr(scale(spatC_db$mass_ratio), "scaled:scale"))+attr(scale(spatC_db$mass_ratio), "scaled:center"),p_disreal= (as.numeric(as.character(group))*attr(scale(spatC_db$p_distance), "scaled:scale"))+attr(scale(spatC_db$p_distance), "scaled:center")) %>%mutate(across(p_disreal, round, 2))SC_pred_massxl_F <-ggemmeans(SC_mod1_F, terms =c("mass_ratio[all]", "Lat_abs[SCp_latvec_F]"),ci_level=.85,ci.lvl = .85) %>%mutate( mass_real= (x*attr(scale(spatC_db$mass_ratio), "scaled:scale"))+attr(scale(spatC_db$mass_ratio), "scaled:center"),lat_real= (as.numeric(as.character(group))*attr(scale(spatC_db$Lat_abs), "scaled:scale"))+attr(scale(spatC_db$Lat_abs), "scaled:center")) %>%mutate(across(lat_real, round, 2))# Temporal allTCp_distvec <-c(max(tempC_modsdf$p_distance), mean(tempC_modsdf$p_distance),min(tempC_modsdf$p_distance))TC_pred_massxp <-ggemmeans(TC_mod3, terms =c("mass_ratio[-1.3:3.1 by=.2]", "p_distance[TCp_distvec]"),ci_level=.85,ci.lvl = .85) %>%mutate( mass_real= (x*attr(scale(tempC_db$mass_ratio), "scaled:scale"))+attr(scale(tempC_db$mass_ratio), "scaled:center"),p_disreal= (as.numeric(as.character(group))*attr(scale(tempC_db$p_distance), "scaled:scale"))+attr(scale(tempC_db$p_distance), "scaled:center")) %>%mutate(across(p_disreal, round, 2))TC_pred_lat <-ggemmeans(TC_mod1, terms =c("Lat_abs[-1.6:2.3 by=.2]"),ci_level=.85,ci.lvl = .85) %>%mutate(lat_real= (x*attr(scale(tempC_db$Lat_abs), "scaled:scale"))+attr(scale(tempC_db$Lat_abs), "scaled:center"))# Temporal FelidaeTC_pred_pxdiet_F1 <-ggemmeans(TC_mod1_F, terms =c("p_distance[-0.16:1.28 by=.2]", "diet_dist[Diff_diet]"),ci_level=.85,ci.lvl = .85) %>%mutate(p_real= (x*attr(scale(tempC_db$p_distance), "scaled:scale"))+attr(scale(tempC_db$p_distance), "scaled:center"))TC_pred_pxdiet_F2 <-ggemmeans(TC_mod1_F, terms =c("p_distance[-2.3:1.08 by=.2]", "diet_dist[Same_diet]"),ci_level=.85,ci.lvl = .85) %>%mutate(p_real= (x*attr(scale(tempC_db$p_distance), "scaled:scale"))+attr(scale(tempC_db$p_distance), "scaled:center"))TC_pred_pxdiet_F <-rbind(TC_pred_pxdiet_F1, TC_pred_pxdiet_F2)# Temporal Canidae TC_pred_mass_C <-ggemmeans(TC_mod1_C, terms =c("mass_ratio[-1.2:1.28 by=.2]"),ci_level=.85,ci.lvl = .85) %>%mutate(mass_real= (x*attr(scale(tempC_db$mass_ratio), "scaled:scale"))+attr(scale(tempC_db$mass_ratio), "scaled:center"))TC_pred_lats_C <-ggemmeans(TC_mod1_C, terms =c("Lat_abs[-1.4:2.3]"),ci_level=.85,ci.lvl = .85) %>%mutate(lat_real= (x*attr(scale(tempC_db$Lat_abs), "scaled:scale"))+attr(scale(tempC_db$Lat_abs), "scaled:center"))# Temporal MustelidaeTC_pred_pdist_M <-ggemmeans(TC_mod1_M, terms =c("p_distance[-1.7:1.3 by=.2]"),ci_level=.85,ci.lvl = .85)%>%mutate(p_real= (x*attr(scale(tempC_db$p_distance), "scaled:scale"))+attr(scale(tempC_db$p_distance), "scaled:center"))TC_pred_diet_M <-ggemmeans(TC_mod1_M, terms =c("diet_dist"),ci_level=.85,ci.lvl = .85)TC_pred_lat_M <-ggemmeans(TC_mod1_M, terms =c("Lat_abs[-1.4:2.3 by=.2]"))%>%mutate(lat_real= (x*attr(scale(tempC_db$Lat_abs), "scaled:scale"))+attr(scale(tempC_db$Lat_abs), "scaled:center"))Lat_TC <-bind_rows(TC_pred_lat, TC_pred_lats_C, TC_pred_lat_M) %>%mutate(data_from=c(rep("All species", length(TC_pred_lat$x)),rep("Canidae", length(TC_pred_lats_C$x)),rep("Mustelidae", length(TC_pred_lat_M$x))))```### Plots```{r}# Spatial all data(SC_mass_plot <-ggplot(SC_pred_mass)+geom_hline(yintercept =1, linetype="dashed", size=1)+geom_ribbon(aes(x= mass_real, y= predicted,ymin= conf.low, ymax= conf.high) , alpha=0.7, fill="#fde725")+geom_line(aes(x= mass_real, y=predicted),size=0.8 )+labs(x="ln(Mass ratio)",y="Spatial Overlap",tag="A")+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto")))(SC_pdist_plot <-ggplot(SC_pred_pdist)+geom_hline(yintercept =1, linetype="dashed", size=1)+geom_ribbon(aes(x= p_real, y= predicted,ymin= conf.low, ymax= conf.high) , alpha=0.7, fill="#fde725")+geom_line(aes(x= p_real, y=predicted),size=0.8 )+labs(x="Phylogenetic distance",y="Spatial Overlap",tag="B")+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto")))# Spatial Felidae(SC_mxp_plot_F <-ggplot(SC_pred_massxp_F)+geom_hline(yintercept =1, linetype="dashed", size=1)+geom_ribbon(aes(x= mass_real, y= predicted,ymin= conf.low, ymax= conf.high,fill=as.factor(p_disreal),group=as.factor(p_disreal)) , alpha=0.7, linewidth=1)+scale_fill_viridis_d()+geom_line(aes(x= mass_real, y=predicted, group=as.factor(p_disreal),linetype=as.factor(p_disreal)),size=0.8 )+labs(x="ln(Mass ratio)",y="Spatial Overlap",group="Phylogenetic \ndistance",linetype="Phylogenetic \ndistance",fill="Phylogenetic \ndistance",tag="C")+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto"),legend.direction="horizontal",legend.position =c(0.5, 0.10),legend.background =element_blank()))(SC_mxl_plot_F <-ggplot(SC_pred_massxl_F)+geom_hline(yintercept =1, linetype="dashed", size=1)+geom_ribbon(aes(x= mass_real, y= predicted,ymin= conf.low, ymax= conf.high,fill=as.factor(lat_real),group=as.factor(lat_real)) , alpha=0.7, linewidth=1)+scale_fill_viridis_d()+geom_line(aes(x= mass_real, y=predicted, group=as.factor(lat_real),linetype=as.factor(lat_real)),size=0.8 )+labs(x="ln(Mass ratio)",y="Spatial Overlap",group="Absolute \nlatitude",linetype="Absolute \nlatitude",fill="Absolute \nlatitude",tag="D")+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto"),legend.direction="horizontal",legend.position =c(0.5, 0.10),legend.background =element_blank()))# Temporal all(TC_mxp_plot <-ggplot(TC_pred_massxp)+geom_ribbon(aes(x= mass_real, y= predicted,ymin= conf.low, ymax= conf.high,fill=as.factor(p_disreal),group=as.factor(p_disreal)) , alpha=0.7, linewidth=1)+scale_fill_viridis_d()+geom_line(aes(x= mass_real, y=predicted, group=as.factor(p_disreal),linetype=as.factor(p_disreal)),size=0.8 )+labs(x="ln(Mass ratio)",y="Temporal Overlap",group="Phylogenetic \ndistance",fill="Phylogenetic \ndistance",linetype="Phylogenetic \ndistance",tag="E")+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto"),legend.direction="horizontal",legend.position =c(0.45, 0.10),legend.background =element_blank()))# TEmporal Felidae(TC_pxd_plot_F <-ggplot(TC_pred_pxdiet_F)+geom_ribbon(aes(x= p_real, y= predicted,ymin= conf.low, ymax= conf.high,fill= group, group= group,linetype= group) , alpha=0.7)+scale_fill_viridis_d()+geom_line(aes(x= p_real, y=predicted, group= group,linetype= group),size=0.8 )+labs(x="Phylogenetic distance",y="Temporal Overlap",group="Diet \ndistance",fill="Diet \ndistance",linetype="Diet \ndistance",tag="F")+scale_x_continuous(expand =c(0,0))+scale_y_continuous(limits =c(0.3,1))+theme_bw()+theme(text =element_text(size=13, family ="Roboto"),legend.direction="horizontal",legend.position =c(0.55, 0.10),legend.background =element_blank()))# Temporal Canidae(TC_mass_plot_C <-ggplot(TC_pred_mass_C)+geom_ribbon(aes(x= mass_real, y= predicted,ymin= conf.low, ymax= conf.high) , alpha=0.7, fill="#fde725")+geom_line(aes(x= mass_real, y=predicted),size=0.8 )+labs(x="ln(Mass ratio)",y="Temporal Overlap",tag="G")+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto")))# Temporal Mustelidae (TC_p_plot_M <-ggplot(TC_pred_pdist_M)+geom_ribbon(aes(x= p_real, y= predicted,ymin= conf.low, ymax= conf.high) , alpha=0.7, fill="#fde725")+geom_line(aes(x= p_real, y=predicted),size=0.8 )+labs(x="Phylogenetic distance",y="Temporal Overlap",tag="H")+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto")) )(TC_diet_plot_M <-ggplot(TC_pred_diet_M)+geom_errorbar(aes(x= x, y= predicted,ymin= conf.low, ymax= conf.high) , size=1, width=0.1)+geom_point(aes(x= x, y=predicted),size=4)+labs(x="Diet distance",y="Temporal Overlap",tag="I")+theme_bw()+theme(text =element_text(size=13, family ="Roboto")) )(TC_Lat_plot <-ggplot(Lat_TC)+geom_ribbon(aes(x= lat_real, y= predicted,ymin= conf.low, ymax= conf.high,group= data_from, fill=data_from) , alpha=0.7)+geom_line(aes(x= lat_real, y=predicted, linetype= data_from, group= data_from),size=0.8)+scale_fill_viridis_d()+labs(x="Absolute latitude",y="Temporal Overlap",fill="Dominant \ncompetitor",linetype="Dominant \ncompetitor",tag="J")+theme_bw()+scale_x_continuous(expand =c(0,0))+theme_bw()+theme(text =element_text(size=13, family ="Roboto"),legend.direction="vertical",legend.position =c(1.2, 0.5),legend.background =element_blank()))``````{r}#| message: false#| warning: false#| fig-width: 11.5#| fig-height: 11.5(CPrediction_plot <- (SC_mass_plot + SC_pdist_plot+ SC_mxp_plot_F + SC_mxl_plot_F+ TC_mxp_plot+ TC_pxd_plot_F+ TC_mass_plot_C+ TC_p_plot_M+ TC_diet_plot_M+ TC_Lat_plot)+plot_layout(ncol =3))``````{r}#| eval: falseggsave(filename ="Figs/Cpreds_plot.svg", plot = CPrediction_plot,width =11.5, height =11.5)```
